Estimating Vaccine Efficacy Over Time After a Randomized Study is Unblinded
aa r X i v : . [ s t a t . A P ] F e b Estimating Vaccine Efficacy Over Time After a Randomized Study isUnblinded
Anastasios A. Tsiatis ∗ and Marie Davidian ∗∗ Department of Statistics, North Carolina State University, Raleigh, NC, USA ∗ [email protected] ∗∗ [email protected] Abstract
The COVID-19 pandemic due to the novel coronavirus SARS CoV-2 has inspired remarkable break-throughs in development of vaccines against the virus and the launch of several phase 3 vaccine trials inSummer 2020 to evaluate vaccine efficacy (VE). Trials of vaccine candidates using mRNA delivery systemsdeveloped by Pfizer-BioNTech and Moderna have shown substantial VEs of 94-95%, leading the US Foodand Drug Administration to issue Emergency Use Authorizations and subsequent widespread adminis-tration of the vaccines. As the trials continue, a key issue is the possibility that VE may wane over time.Ethical considerations dictate that all trial participants be unblinded and those randomized to placebobe offered vaccine, leading to trial protocol amendments specifying unblinding strategies. Crossover ofplacebo subjects to vaccine complicates inference on waning of VE. We focus on the particular featuresof the Moderna trial and propose a statistical framework based on a potential outcomes formulationwithin which we develop methods for inference on whether or not VE wanes over time and estimation ofVE at any post-vaccination time. The framework clarifies assumptions made regarding individual- andpopulation-level phenomena and acknowledges the possibility that subjects who are more or less likely tobecome infected may be crossed over to vaccine differentially over time. The principles of the frameworkcan be adapted straightforwardly to other trials.
Key words:
Crossover; Inverse probability weighting; Potential outcomes; Randomized phase 3 vaccinetrial; Waning vaccine efficacy 1
Introduction
The primary objective of a vaccine trial is to estimate vaccine efficacy (VE). Typically, these trials aredouble-blind, placebo-controlled studies in which participants are randomized to either vaccine or placeboand followed for the primary endpoint, which is often time to viral infection, on which inference on VEis based, where VE is defined as a measure of reduction in infection risk for vaccine relative to placebo,expressed as a percentage.Vaccine trials have become the focus of immense global interest as a result of the COVID-19 diseasepandemic due to the novel coronavirus SARS-CoV-2. The pandemic inspired unprecedented scientificbreakthroughs in the rapid development of vaccines against SARS-CoV-2, culminating in the launch ofseveral large phase 3 vaccine trials in Summer 2020. Trials in the US studying the vaccine candidatesusing messenger RNA (mRNA) delivery systems developed by Pfizer-BioNTech and Moderna began inJuly 2020 and demonstrated substantial evidence of VEs of 94-95% at interim analyses, leading the USFood and Drug Administration (FDA) to issue Emergency Use Authorizations (EUAs) for both vaccinesin December 2020 and to the rollout of vaccination programs shortly thereafter.Implicit in the primary analysis in these trials is the assumption that VE is constant over the studyperiod and, with primary endpoint time to infection, VE is represented by the 1 − the ratio of the hazardrate for vaccine to that for placebo, estimated based on a Cox proportional hazards model. As the trialscontinue following the EUAs, among the many issues to be addressed is the possibility that VE may waneover time. Principled evaluation of the nature and extent of waning of VE is of critical public healthimportance, as waning has implications for measures to control the pandemic. Were all participants in thetrials to continue on their randomized assignments (vaccine or placebo), evaluation of potential waningof VE would be straightforward. However, once efficacy is established, ethical considerations dictate thepossibility of unblinding all participants and offering the vaccine to those randomized to placebo. Afterconsultation with stakeholders, Pfizer and Moderna issued amendments to their trial protocols specifyingunblinding strategies and modifications to planned analyses.Crossover of placebo subjects to vaccine of necessity complicates inference on waning of VE andhas inspired recent research (Follmann et al., 2020; Fintzi and Follmann, 2021; Lin, Zeng and Gilbert,2021). We propose a statistical framework within which we develop methods for inference on whetheror not VE wanes over time based on data where subjects are unblinded and those on placebo may crossover to vaccine and in which assumptions made regarding individual and population phenomena aremade transparent. It is possible that subjects who are more or less likely to become infected could beunblinded and cross over to vaccine differentially over time, which could lead to biased inferences due to2onfounding; accordingly, this possibility is addressed explicitly in the framework. The first author (AAT)has the privilege of serving on the Data and Safety Monitoring Board for all US government-sponsoredCOVID-19 vaccine trials and is thus well-acquainted with the unblinding approach for the Moderna trial.Accordingly, the development is based on the specifics of this trial, but the principles can be adapted tothe features of other trials.In Section 2, we review the Moderna trial and the resulting data. We present a conceptual frameworkin which we precisely define VE as a function of time post-vaccination in Section 3. In Section 4, wedevelop a formal statistical framework within which we propose methodology for estimation of VE anddescribe its practical implementation in Section 5. Simulations demonstrating performance are presentedin Section 6. We first describe the timeline of the Moderna Coronavirus Efficacy (COVE) trial (Baden et al., 2020) onthe scale of calendar time. The trial opened on July 27, 2020 (time 0), and reached full accrual at time T A (October 23, 2020). On December 11, 2020, denoted T P , the FDA issued an EUA for the Pfizer vaccine,followed by an EUA for the Moderna mRNA-1273 vaccine on T M = December 18, 2020. Amendment 6of the study protocol was issued on December 23, 2020 and specified the unblinding strategy (see Figure2 of the protocol) under which, starting on T U = December 24, 2020, study participants are scheduledon a rolling basis over several months for Participant Decision clinic visits (PDCVs) at which they willbe unblinded. If originally randomized to vaccine, participants continue to be followed; if randomizedto placebo, participants can receive the Moderna vaccine or refuse. Let T C denote the time at whichall PDCVs have taken place. The study will continue until time T F at which all participants will havecompleted full follow-up at 24 months after initial treatment assignment. Assume that the analysis ofvaccine efficacy using the methods in Sections 4.4 and 5 takes place at time T C ≤ L ≤ T F , where allparticipants have achieved the primary endpoint, requested to be unblinded, or attended the PDCV by L . Under this scheme, we characterize the data on a given participant as follows. Let 0 ≤ E ≤ T A denote the calendar time at which the subject entered the trial, X denote baseline covariates, and A = 0(1) if assigned to placebo (vaccine). Denote observed time to infection on the scale of calendar time as U , and ∆ = I( U ≤ L ), where I( B ) = 1 if B is true and 0 otherwise. At T P , availability of the Pfizervaccine commenced, at which point some subjects not yet infected requested to be unblinded. Denoteby R (calendar time) the minimum of (i) time to such an unblinding, in which case T P ≤ R < T U , and3 able 1: Summary of notation. All times are on the scale of calendar time, where time 0 is the start of the trial. Trial Milestones T A Full accrual reached, October 23, 2020 T P Pfizer granted EUA, December 11, 2020 T M Moderna granted EUA, December 18, 2020 T U Participant Decision clinic visits (PDCVs) commence, December 24, 2020 T C PDCVs conclude T F Follow-up concludes, trial ends ℓ Lag between initial vaccine dose and full efficacy, 6 weeks, T P − T A > ℓL Time of analysis of vaccine efficacy using the proposed methods;
L > time atwhich all subjects have achieved the endpoint, requested unblinding, or attendedthe PDCV, L ≤ T F Observed Data on a Trial Participant E Study entry time, 0 ≤ E ≤ T A X Baseline information A Treatment assignment, placebo, A = 0, or vaccine, A = 1 U, ∆ Time to infection, indicator of infection by time L , ∆ = I( U ≤ L ) R, Γ Time to requested unblinding, PDCV/requested unblinding, or infection, whichevercomes firstΓ = 0: R = U , infection occurs before requested/offered unblindingΓ = 1: R = time to requested unblinding, T P ≤ R < T U Γ = 2: R = time to PDCV or requested unblinding, T U ≤ R < T C Ψ If A = 0, Γ ≥
1, indicator or whether subject receives Moderna vaccine, Ψ = 1,or refuses, Ψ = 0define Γ = 1; (ii) time of PDCV, so T U ≤ R < T C , and let Γ = 2; or (iii) time to infection, in which case R = U and Γ = 0. If Γ ≥ A = 1, so that the subject was randomized to vaccine, s/he continuesto be followed; if A = 0, s/he can choose to receive the Moderna vaccine, Ψ = 1 or refuse, Ψ = 0. Wedistinguish the cases Γ = 1 and 2 to acknowledge different unblinding dynamics before and after T U .Because a very small number of participants requested unblinding before T P , and, although the protocolallows participants to refuse unblinding at PDCV, all subjects are strongly encouraged to unblind, we donot include these possibilities in the formulation.Table 1 summarizes the timeline and observed data. The trial data are thus O i = { E i , X i , A i , U i , ∆ i , R i , Γ i , I(Γ i ≥ , A i = 0)Ψ i } , i = 1 , . . . , n, (1)independent and identically distributed (iid) across i .4 Conceptualization of Vaccine Efficacy
Similar to Halloran, Longini, and Struchiner (1996) and Longini and Halloran (1996), we consider thefollowing framework in which to conceptualize vaccine efficacy. The study population, comprising indi-viduals for which inference on vaccine efficacy is of interest, is that of individuals susceptible to infection,represented by the trial participants. There is a population of individuals outside the trial with which trialparticipants interact, assumed to be much larger than the number of participants, so that interactionsamong participants are much less likely than interactions with the outside population. The probabilitythat a trial participant will become infected at calendar time t depends on three factors: c ( t ), the contactrate, the number of contacts with the outside population per unit time; p ( t ), the prevalence of infectionsin the outside population at t ; and π ( t ), the transmission probability at t , the probability a susceptibleindividual in the study population will become infected per contact with an infected individual from theoutside population. Dependence of π ( t ) on time acknowledges the emergence of new variants of the virus,which may be be more or less virulent, as in the COVID-19 pandemic. Assuming random mixing, p ( t ) c ( t )is the contact rate at time t with infected individuals, and the infection rate at time t is p ( t ) c ( t ) π ( t ).We adapt this framework to the COVID-19 pandemic. The prevalence rate in the pandemic can varysubstantially in time and space, so denote by S the trial site at which a participant is enrolled, andlet p ( t, s ) be the prevalence at time t at site S = s . Although p ( t, s ) varies by t and s , assume it isunaffected by the individuals in the trial and thus represents an external force. We view the contact rateas individual specific; accordingly, for an arbitrary individual in the study population, let the randomvariables { c b ( t ) , c b ( t ) , c u ( t ) , c u ℓ ( t ) , c u ( t ) } denote potential contact rates. These potential outcomes can beregarded as individual-specific behavioral characteristics of trial participants, where some may be morecareful and make fewer contacts while others take more risks, and behavior can vary over time and byvaccination and blinding status. Here, c ba ( t ) is the contact rate at time t if the individual were to receivevaccine, a = 1, or placebo, a = 0, and be blinded to this assignment; by virtue of blinding, it is reasonableto take c b ( t ) = c b ( t ) = c b ( t ). The Moderna vaccine is administered in two doses, ideally 4 weeks apart, andis not thought to achieve full efficacy until 2 weeks following the second dose. Thus, letting ℓ denote thelag between initial dose and full efficacy, c u ℓ ( t ) and c u ( t ) reflect behavior of an individual who is unblindedand vaccinated in the periods prior to ℓ and after ℓ , respectively, allowing for unblinded vaccinees to, e.g.,behave more cautiously before full efficacy is achieved. The rate c u ( t ) reflects behavior of an unblindedindividual on placebo and does not play a role in the development. Similar to the stable unit treatmentvalue assumption (Rubin, 1980), assume that c u ℓ ( t ) and c u ( t ) are the same whether the individual wasrandomized to vaccine and unblinded before t or was randomized to placebo and subsequently unblinded5nd crossed over to vaccine before t .Finally, for an arbitrary participant, let the random variable π ( t ) be the potential individual-specifictransmission probability per contact at t if s/he were to receive placebo, and let π ( t, τ ) be the same ifs/he were to receive vaccine and have been vaccinated for τ ≥ τ and thus consider whether or not VE wanesover time since vaccination.With the set of potential outcomes for an arbitrary individual in the study population who enrolls atsite S thus given by { c b ( t ) , c u ( t ) , c u ℓ ( t ) , c u ( t ) t > , π ( t ) , π ( t, τ ) , τ ≥ } , the infection rate in the studypopulation at calendar time t if all individuals were to receive placebo and be blinded to that assignmentis I b ( t ) = E { p ( t, S ) c b ( t ) π ( t ) } ; likewise, the infection rate at t if all individuals were to receive vaccine attime t − τ and be blinded to that assignment is I b ( t, τ ) = E { p ( t, S ) c b ( t ) π ( t, τ ) } . The relative infectionrate at t is then R b ( t, τ ) = I b ( t, τ ) I b ( t ) = E { p ( t, S ) c b ( t ) π ( t, τ ) } E { p ( t, S ) c b ( t ) π ( t ) } . (2)Accordingly, vaccine efficacy at time t after vaccination at t − τ is V E ( t, τ ) = 1 − R b ( t, τ ), reflectingthe proportion of infections at t that would be prevented if the study population were vaccinated and onvaccine for τ units of time during the blinded phase of the study.In the sequel, we assume that R b ( t, τ ) and thus V E ( t, τ ) depend only on τ and write R b ( τ ) and V E ( τ ) = 1 − R b ( τ ). This assumption embodies the belief that, although infection rates may changeover time, the relative effect of vaccine to placebo remains approximately constant and holds if (i) { π ( t, τ ) , π ( t ) } ⊥⊥ { S, c b ( t ) }| X , where ⊥⊥ means “independent of” and this independence is conditionalon X ; and (ii) E { π ( t, τ ) | X } /E { π ( t ) | X } = q ( τ ), so does not depend on t and X . Condition (i) reflectsthe interpretation of π ( t, τ ) and π ( t ) as inherent biological characteristics of an individual, whereas S and c b ( t ) are external and behavioral characteristics, respectively; thus, once common individual andexternal baseline covariates are taken into account, biological and geographic/behavioral characteristicsare unrelated. Condition (ii) implies that, although new viral variants may change transmission prob-abilities under both vaccine and placebo over time, this change stays in constant proportion, and thisproportion is similar for individuals with different characteristics. Further discussion is given in Section 7and Appendix B.Within this framework, the goal of inference on waning of VE based on the data from the trial can bestated precisely as inference on V E ( τ ) = 1 − R b ( τ ), τ ≥ ℓ , so reflecting VE after full efficacy is achieved.It is critical to recognize that, like estimands of interest in most clinical trials, V E ( τ ) represents VEat time since vaccination τ under the original conditions of the trial, under which all participants are6linded. The challenge we address in subsequent sections is how to achieve valid inference on V E ( τ ), τ ≥ ℓ , using data from the modified trial in which blinded participants are unblinded in a staggeredfashion, with placebo subjects offered the option to receive vaccine.We propose a semiparametric model within which we cast this objective. Let I u ℓ ( t, τ ) = E { p ( t, S ) c u ℓ ( t ) π ( t, τ ) } , τ < ℓ , and I u ( t, τ ) = E { p ( t, S ) c u ( t ) π ( t, τ ) } , τ ≥ ℓ , be the infection rates in the study population at t if all individuals were to receive vaccine at time t − τ and be unblinded to that fact. Analogous to (i)above, assume that { π ( t, τ ) , π ( t ) } ⊥⊥ { S, c u ℓ ( t ) , c u ( t ) }| X , and continue to assume condition (ii). Then,for two values τ , τ of τ , it is straightforward that (see Appendix A) I u ℓ ( t, τ ) I u ℓ ( t, τ ) = R b ( τ ) R b ( τ ) , τ , τ < ℓ ; I u ( t, τ ) I u ( t, τ ) = R b ( τ ) R b ( τ ) , τ , τ ≥ ℓ. (3)Defining I u ℓ ( t ) = I u ℓ ( t,
0) = E { p ( t, S ) c u ℓ ( t ) π ( t, } and I u ( t ) = I u ( t, ℓ ) = E { p ( t, S ) c u ( t ) π ( t, ℓ ) } , by (3)with τ = τ and τ = 0 ( ℓ ) on the left (right) hand side, the infection rates at t if all individuals in thestudy population were unblinded and to receive vaccine at time t − τ are I u ℓ ( t, τ ) = I u ℓ ( t ) R b ( τ ) R b (0) , τ < ℓ ; I u ( t, τ ) = I u ( t ) R b ( τ ) R b ( ℓ ) , τ ≥ ℓ. (4)Likewise, from (2), the infection rate at t if all individuals in the study population were blinded and toreceive vaccine at time t − τ is I b ( t, τ ) = I b ( t ) R b ( τ ) . (5)We now represent the infection rate ratio R b ( τ ) as R b ( τ ; θ ) = exp { ζ ( τ ) } I( τ < ℓ ) + exp { θ + g ( τ − ℓ ; θ ) } I( τ ≥ ℓ ) , θ = ( θ , θ T ) T , (6)where ζ ( τ ) is an unspecified function of τ ; θ and θ are real- and vector-valued parameters, respectively;and g ( u ; θ ) is a real-valued function of such that g (0; θ ) = 0 for all θ and g ( u ; 0) = 0. For example,taking g ( u ; θ ) = θ u yields R b ( τ ; θ ) = exp { θ + θ ( τ − ℓ ) } , τ ≥ ℓ , in which case θ = 0 implies that V E ( τ ) = 1 − R b ( τ ), τ ≥ ℓ , does not change with time since vaccination, and θ > V E ( τ )decreases with increasing τ ; i.e., exhibits waning. More complex specifications of g ( u ; θ ) using splines(e.g., Fintzi and Follmann, 2021) or piecewise constant functions could be made; e.g., for v < v ≤ L , g ( u ; θ ) = θ I( v < u ≤ v ) + θ I( u > v ) , θ = ( θ , θ ) T . (7)Under this model, (5) and (4) can be written as I b ( t, τ ) = I b ( t ) (cid:2) exp { ζ ( τ ) } I( τ < ℓ ) + exp { θ + g ( τ − ℓ ; θ ) } I( τ ≥ ℓ ) (cid:3) , I u ℓ ( t, τ ) = I u ℓ ( t ) exp { ζ ( τ ) } , τ < ℓ, I u ( t, τ ) = I u ( t ) exp { g ( τ − ℓ ; θ ) } , τ ≥ ℓ. (8)Thus, to estimate V E ( τ ) for any τ and make inference on potential waning of VE, we must develop aprincipled approach to estimation of θ based on the data from the modified trial in which participantsare unblinded and those on placebo may cross over to vaccine.7 Statistical Framework
Estimation of
V E ( τ ), equivalently R b ( τ ), would be straightforward for any τ ≥ ℓ over the entire follow-upperiod if all participants remained on their assigned treatments throughout the trial. However, subjectsrandomized to placebo have the option to cross over to vaccine on or after T P . For τ < T P , it is possibleto estimate R b ( τ ) because, due to randomization, for t < T P we have representative samples of blindedsubjects on vaccine and placebo and thus information on I b ( t, τ ) and I b ( t ), so can estimate θ andcomponents of θ identified for such τ ; e.g., in (7) depending on the values of v and v . At T P ≤ t < T C ,the data comprise a mixture of blinded and unblinded participants, where, within the latter group, thoseon placebo may have crossed over to vaccine. Here, information, albeit diminishing during the interval[ T P , T C ), on I b ( t, τ ) and I b ( t ) is available from those participants not yet unblinded, which contributesto estimation of θ and components of θ . Information is also available on I u ( t, τ ) from individuals whowere originally randomized to vaccine and provide information on longer τ , and from individuals whorecently crossed over to vaccine and provide information on shorter τ . For t ≥ T C , there are no longerblinded participants, so that information is available only on I u ( t, τ ). For these latter groups, for longer τ ≥ ℓ and shorter τ ≥ ℓ , I u ( t, τ ) / I u ( t, τ ) = exp (cid:2) g { τ − ℓ ; θ } − g { τ − ℓ ; θ } (cid:3) , and, because of themixture of times since vaccination, θ can be fully estimated.Through the following potential outcomes formulation and under suitable assumptions, in the nextseveral sections we develop an approach to estimation of θ based on the observed data (1) that embodiesthe foregoing intuitive principles. Denote by T ∗ ( e, r ) the potential time to infection on the scale of patient time for an arbitrary individualin the study population if s/he were to enter the trial at calendar time e , receive placebo and be blindedto that fact, and, if not infected by calendar time r , be unblinded and cross over to vaccine at r . Let T ∗ ( e ) = T ∗ ( e, ∞ ), if s/he is never crossed over to receive vaccine. Similarly, define T ∗ ( e, r ) to be thepotential time to infection (patient time scale) for an arbitrary individual if s/he were to enter the trialat e , receive vaccine and be blinded to that fact, and, if not infected by r , be unblinded at r ; anddefine T ∗ ( e ) = T ∗ ( e, ∞ ). We make the consistency assumptions that T ∗ ( e, r ) = T ∗ ( e ) if T ∗ ( e ) < r and T ∗ ( e, r ) = T ∗ ( e ) if T ∗ ( e ) < r . For a = 0 ,
1, denote the hazard at calendar time t , t > e , by λ a ( t, e, r ) = lim dt → pr { t ≤ T ∗ a ( e, r ) + e < t + dt | T ∗ a ( e, r ) + e ≥ t } , a = 0 , , (9)8here the addition of e induces a shift from patient to calendar time. Denote the set of all potentialoutcomes as W ∗ = { T ∗ ( e, r ) , T ∗ ( e, r ); e > , r > e } . The development in Section 3 is in terms of infection rates at the individual-specific and populationlevels. Population-level hazard rates such as (9) are not equivalent to population-level infection rates.However, we argue in Appendix C that, because the probabilities of infection under vaccine and placeboduring the course of the trial are small, population-level hazard rates and population-level infectionrates are approximately equivalent; this assumption is implicit in the standard primary analysis notedin Section 1. Thus, to reflect this, we use familiar notation and write λ b ( t ) = I b ( t ), λ uℓ ( t ) = I u ℓ ( t ), and λ u ( t ) = I u ( t ). Under these conditions, using (8), we can write for t > eλ ( t, e, r ) = λ b ( t )I( t < r ) + λ uℓ ( t ) exp { ζ ( t − r ) } I( t − r < ℓ )+ λ u ( t ) exp { g ( t − r − ℓ ; θ ) } I( t − r ≥ ℓ ) , (10) λ ( t, e, r ) = λ b ( t ) (cid:2) exp { ζ ( t − e ) } I( t − e < ℓ ) + exp { θ + g ( t − e − ℓ ; θ ) } I( t − e ≥ ℓ ) (cid:3) I( t < r )+ λ u ( t ) exp { g ( t − e − ℓ ; θ ) } I( t ≥ r ) , (11)where (11) follows because r ≥ T P , e ≤ T A , T P − T A > ℓ . Define the counting processes for infectionby N ∗ a ( t, e, r ) = I { T ∗ a ( e, r ) + e ≤ t } and N ∗ a ( t, e ) = N ∗ a ( t, e, ∞ ), and the at-risk processes by Y ∗ a ( t, e, r ) =I { T ∗ a ( e, r ) + e ≥ t } and Y ∗ a ( t, e ) = Y ∗ a ( t, e, ∞ ), a = 0 , t < r , then N ∗ a ( t, e, r ) = N ∗ a ( t, e ), Y ∗ a ( t, e, r ) = Y ∗ a ( t, e ), a = 0 ,
1. For a = 0 ,
1, let Λ a ( t, e, r ) = R t λ a ( u, e, r ) du be the cumulative hazard. Because E { dN ∗ a ( t, e, r ) | Y ∗ a ( t, e, r ) } = d Λ a ( t, e, r ) Y ∗ a ( t, e, r ), a = 0 ,
1, it follows that { dN ∗ a ( t, e, r ) − d Λ a ( t, e, r ) Y ∗ a ( t, e, r ) } , a = 0 ,
1, are mean-zerocounting process increments. Thus, any linear combination of these increments over t, e, r can be usedto define unbiased estimating functions in W ∗ of quantities of interest. In Appendix D, we formulate aparticular set of estimating functions such that, given iid potential outcomes W ∗ i , i = 1 , . . . , n , lead toconsistent and asymptotically normal estimators for { Λ b ( t ) , Λ u ( t ) , θ T } T , Λ k ( t ) = R t λ k ( u ) du , k = b, u .Because interest focuses on V E ( τ ) for τ ≥ ℓ , estimation of Λ uℓ ( t ) = R t λ uℓ ( u ) du and ζ ( · ) is not considered.For fixed t , 0 ≤ t ≤ L , the estimating functions for Λ b ( t ) and Λ u ( t ) are, respectively, E ∗ Λ b { W ∗ ; Λ b ( t ) , θ } = I( t < T C ) Z min( t, T A )0 { dN ∗ ( t, e ) − d Λ b ( t ) Y ∗ ( t, e ) } e w ( t, e ) de (12)+I( t ≥ ℓ ) Z min( t − ℓ, T A )0 (cid:2) dN ∗ ( t, e ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ )I( t − e ≥ ℓ ) } Y ∗ ( t, e ) (cid:3) e w ( t, e ) de ! , ∗ Λ u { W ∗ ;Λ u ( t ) , θ } = I( t ≥ T P + ℓ ) Z T A Z min( t − ℓ, T C ) T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − r − ℓ ; θ )I( t − r ≥ ℓ ) } Y ∗ ( t, e, r ) (cid:3) w ( t, e, r ) dr de ! + I( t ≥ T P ) Z T A Z min( t, T C ) T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − e − ℓ ; θ ) }× Y ∗ ( t, e, r ) (cid:3) w ( t, e, r )I( t ≥ r ) dr de ! , (13)where e w a ( t, e ) and w a ( t, e, r ), a = 0 ,
1, are arbitrary nonnegative weight functions, specification of whichis discussed later. The estimating function for θ is given by E ∗ θ { W ∗ ; Λ b ( · ) , Λ u ( · ) , θ } = Z T C ℓ Z min( t − ℓ, T A )0 g θ ( t − e − ℓ ) (cid:2) dN ∗ ( t, e ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ ) × I( t − e ≥ ℓ ) } Y ∗ ( t, e ) (cid:3) e w ( t, e ) de + Z L T P + ℓ Z T A Z min( t − ℓ, T C ) T P g θ ( t − r − ℓ ) (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − r − ℓ ; θ )I( t − r ≥ ℓ ) }× Y ∗ ( t, e, r ) (cid:3) w ( t, e, r ) dr de (14)+ Z L T P Z T A Z min( t, T C ) T P g θ ( t − e − ℓ ) (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − e − ℓ ; θ ) }× Y ∗ ( t, e, r ) (cid:3) w ( t, e, r )I( t ≥ r ) dr de, where g θ ( u ) = ∂/∂θ { g ( u ; θ ) } . Analogous to Yang, Tsiatis, and Blazing (2018), envisioning (12)-(14) ascharacterizing a system of estimating functions E ∗ { W ∗ ; Λ b ( · ) , Λ u ( · ) , θ } = [ E ∗ Λ b { W ∗ ; Λ b ( t ) , θ } , E ∗ Λ u { W ∗ ; Λ u ( t ) , θ } , ≤ t ≤ L, E ∗ θ { W ∗ ; Λ b ( · ) , Λ u ( · ) , θ } T ] T , if we could observe W ∗ i , i = 1 . . . , n , we would estimate d Λ b ( · ) , d Λ u ( · ) , θ by solving the estimating equa-tions P ni =1 E ∗ ( W ∗ i ; Λ b ( · ) , Λ u ( · ) , θ ) } = 0. Of course, the potential outcomes W ∗ i , i = 1 , . . . , n , are not observed. However, we now present assump-tions under which we can exploit the developments in the last section to derive estimating equationsyielding estimators based on the observed data (1).Define the indicator that a participant is observed to be infected at time t by dN ( t ) = I( U = t, ∆ = 1),the observed at-risk indicator at t by Y ( t ) = I( E < t ≤ U ), and I ( t, e ) = (1 − A )I( E = e )I( R ≥ t ) , I ( t, e ) = A I( E = e )I( R ≥ t ) ,I ( t, e, r ) = (1 − A )I( E = e ) { I( R = r, Γ = 1 , Ψ = 1) + I( R = r, Γ = 2 , Ψ = 1) } , (15) I ( t, e, r ) = A I( E = e ) { I( R = r, Γ = 1) + I( R = r, Γ = 2) } . a ( t, e ) = 1 indicates that a subject entering the trial at time e and randomized to placebo ( a = 0) orvaccine ( a = 1) has not yet been infected or unblinded by t . For t > r , I ( t, e, r ) = 1 indicates that asubject randomized to placebo at entry time e is unblinded (either by request or at a PDCV) at time r and crosses over to vaccine at r , and I ( t, e, r ) = 1 if a subject randomized to vaccine at entry time e isunblinded at r . Make the consistency assumptions I a ( t, e ) dN ( t ) = I a ( t, e ) dN ∗ a ( t, e ) , I a ( t, e ) Y ( t ) = I a ( t, e ) Y ∗ a ( t, e ) , a = 0 , ,I ( t, e, r ) dN ( t ) = I ( t, e, r ) dN ∗ ( t, e, r ) , I ( t, e, r ) Y ( t ) = I ( t, e, r ) Y ∗ ( t, e, r ) , (16) I ( t, e, r ) dN ( t ) = I ( t, e, r ) dN ∗ ( t, e, r ) , I ( t, e, r ) Y ( t ) = I ( t, e, r ) Y ∗ ( t, e, r ) . We now make assumptions similar in spirit to those adopted in observational studies. By randomiza-tion, A ⊥⊥ ( X, E, W ∗ ) , (17)where we subsume the site indicator S in X , and let p A = pr( A = 1). It is realistic to assume that the mixof baseline covariates changes over the accrual period; e.g., during the trial, because of lagging accrual ofelderly subjects and subjects from underrepresented groups, an effort was made to increase participationof these groups in the latter part of the accrual period. Accordingly, we allow the distribution of entrytime E to depend on X , and denote its conditional density as f E | X ( e | x ). We make the no unmeasuredconfounders assumption E ⊥⊥ W ∗ | X. (18)Define the hazard functions of unblinding in the periods between the Pfizer EUA and the start ofPDCVs and after the start of PDCVs, respectively, as λ R, ( r | X, A, E, W ∗ ) = lim dr → pr( r ≤ R < r + dr, Γ = 1 | R ≥ r, X, A, E, W ∗ ) , T P ≤ r < T U λ R, ( r | X, A, E, W ∗ ) = lim dr → pr( r ≤ R < r + dr, Γ = 2 | R ≥ r, X, A, E, W ∗ ) , T U ≤ r < T C , where λ R,j ( r | X, A, E, W ∗ ) = 0 for r ≥ T U ( j = 1) and r ≥ T C ( j = 2). Because the accrual periodwas short relative to the length of follow-up, we take these unblinding hazard functions to not dependon E , although including such dependence is straightforward; and, similar to a noninformative censoringassumption, to not depend on W ∗ and write λ R,j ( r | X, A, E, W ∗ ) = λ R,j ( r | X, A ) , j = 1 , . (19)Define K R ( r | X, A ) = exp[ −{ Λ R, ( r | X, A ) + Λ R, ( r | X, A ) } ], Λ R,j ( r | X, A ) = R r T j λ R,j ( u | X, A ) du , T j = T P ( j = 1) or T j = T U ( j = 2). Because λ R, ( r | X, A ) and λ R, ( r | X, A ) are defined on the nonoverlapping11ntervals [ T P , T U ) and [ T U , T C ), respectively, with K R,j ( r | X, A ) = exp {− Λ R,j ( r | X, A ) } , j = 1 , K R ( r | X, A ) = 1 , r < T P , = K R, ( r | X, A ) , T P ≤ r < T U , = K R, ( T U | X, A ) K R, ( r | X, A ) , T U ≤ r < T C , = 0 , r ≥ T C . Finally, define f R,j ( r | X, A ) = K R ( r | X, A ) λ R,j ( r | X, A ), j = 1 , | X, E, Γ , R, W ∗ ) be the probability that a placebo participant unblinded at R agreesto receive the Moderna vaccine. Similar to (19), we assume this probability does not depend on E, W ∗ ;moreover, because the unblinding interval [ T P , T C ) is very short relative to the length of follow-up, weassume it does not depend on R but does depend on the unblinding dynamics at R . Thus, writepr(Ψ = 1 | X, E, Γ , R, W ∗ ) = pr(Ψ = 1 | X, Γ) = p Ψ ( X, Γ) . (20) We now outline, under the assumptions (16)-(20), which we take to hold henceforth, how we can developunbiased estimating equations based on the observed data yielding consistent and asymptotically normalestimators for d Λ b ( · ) , d Λ u ( · ) , θ . The basic premise is to use inverse probability weighting (IPW) toprobabilistically represent potential outcomes in terms of the observed data to mimic the estimatingfunctions (12)-(14).Considering (15), define the inverse probability weights h ( t, e | X ) = (1 − p A ) f E | X ( e | X ) K R ( t | X, A = 0) , h ( t, e | X ) = p A f E | X ( e | X ) K R ( t | X, A = 1) ,h ( e, r | X ) = (1 − p A ) f E | X ( e | X ) × { f R, ( r | X, A = 0) p Ψ ( X, Γ = 1) + f R, ( r | X, A = 0) p Ψ ( X, Γ = 2) } ,h ( e, r | X ) = p A f E | X ( e | X ) { f R, ( r | X, A = 1) + f R, ( r | X, A = 1) } . We show in Appendix E that E (cid:26) I ( t, e ) dN ( t ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e ) , E (cid:26) I ( t, e ) Y ( t ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = Y ∗ ( t, e ) (21) E (cid:26) I ( t, e ) dN ( t ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e ) , E (cid:26) I ( t, e ) Y ( t ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = Y ∗ ( t, e ) , (22) E (cid:26) I ( t, e, r ) dN ( t ) h ( e, r | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e, r ) , E (cid:26) I ( t, e, r ) Y ( t ) h ( e, r | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = Y ∗ ( t, e, r ) , (23) E (cid:26) I ( t, e, r ) dN ( t ) h ( e, r | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e, r ) , E (cid:26) I ( t, e, r ) Y ( t ) h ( e, r | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = Y ∗ ( t, e, r ) . (24)12o obtain observed data analogs to the estimating functions (12)-(14), based on the equalities in (21)-(24), we substitute the IPW expressions in the conditional expectations on the left hand sides. Using(15) and (21)-(22), the analog to (12) is given by E Λ b { O ; Λ b ( t ) , θ } = I( t < T C ) Z min( t, T A )0 I ( t, e ) h ( t, e | X ) { dN ( t ) − d Λ b ( t ) Y ( t ) } e w ( t, e ) de + I( t ≥ ℓ ) Z min( t − ℓ, T A )0 I ( t, e ) h ( t, e | X ) (cid:2) dN ( t ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ ) × I( t − e ≥ ℓ ) } Y ( t ) (cid:3) e w ( t, e ) de ! = I( t < T C ) (cid:18) (1 − A )I( R ≥ t ) h ( t, E | X ) { dN ( t ) − d Λ b ( t ) Y ( t ) } e w ( t, E )+ A I( E + ℓ ≤ t ≤ R ) h ( t, E | X ) (cid:2) dN ( t ) − d Λ b ( t ) exp { θ + g ( t − E − ℓ ; θ ) } Y ( t ) (cid:3) e w ( t, E ) (cid:19) . (25)Likewise, using (23)-(24), the analog to (13) is E Λ u { O ;Λ u ( t ) , θ } = I( t ≥ T P + ℓ ) Z T A Z min( t − ℓ, T C ) T P I ( t, e, r ) h ( e, r | X ) (cid:2) dN ( t ) − d Λ u ( t ) exp { g ( t − r − ℓ ; θ )I( t − r ≥ ℓ ) } Y ( t ) (cid:3) w ( t, e, r ) dr de ! + I( t ≥ T P ) Z T A Z min( t, T C ) T P I ( t, e, r ) h ( e, r | X ) (cid:2) dN ( t ) − d Λ u ( t ) exp { g ( t − e − ℓ ; θ ) }× Y ( t ) (cid:3) w ( t, e, r )I( t ≥ r ) dr de ! = I( t ≥T P + ℓ ) (cid:18) (1 − A )I( t − R ≥ ℓ ) { I(Γ = 1 , Ψ = 1) + I(Γ = 2 , Ψ = 1) } h ( E, R | X ) × (cid:2) dN ( t ) − d Λ u ( t ) exp { g ( t − R − ℓ ; θ ) } Y ( t ) (cid:3) w ( t, E, R ) (cid:19) + I( t ≥ T P ) (cid:18) A I( t > R ) { I(Γ = 1) + I(Γ = 2) } h ( E, R | X ) × (cid:2) dN ( t ) − d Λ u ( t ) exp { g ( t − E − ℓ ; θ ) } Y ( t ) (cid:3) w ( t, E, R ) (cid:19) . (26)A entirely similar representation E θ { O ; Λ b ( · )Λ u ( · ) , θ } of (14) in terms of the observed data can be deducedand is suppressed for brevity. 13o simplify notation, based on (25), (26), and the analogous expression for (14), define d e N b ( t ) = dN ( t ) (cid:26) (1 − A )I( R ≥ t ) e w ( t, E ) h ( t, E | X ) + A I( E + ℓ ≤ t ≤ R ) e w ( t, E ) h ( t, E | X ) (cid:27)e Y b ( t ) = Y ( t ) (cid:20) (1 − A )I( R ≥ t ) e w ( t, E ) h ( t, E | X ) + A I( E + ℓ ≤ t ≤ R ) e w ( t, E ) h ( t, E | X ) exp { θ + g ( t − E − ℓ ; θ ) } (cid:21) d e N u ( t ) = dN ( t ) (cid:20) (1 − A )I( t − R ≥ ℓ ) { I(Γ = 1 , Ψ = 1) + I(Γ = 2 , Ψ = 1) } w ( t, E, R ) h ( E, R | X )+ A I( t > R ) { I(Γ = 1) + I(Γ = 2) } w ( t, E, R ) h ( E, R | X ) (cid:21)e Y u ( t ) = Y ( t ) (cid:20) (1 − A )I( t − R ≥ ℓ ) { I(Γ = 1 , Ψ = 1) + I(Γ = 2 , Ψ = 1) } w ( t, E, R ) h ( E, R | X ) × exp { g ( t − R − ℓ ; θ )+ A I( t > R ) { I(Γ = 1) + I(Γ = 2) } w ( t, E, R ) h ( E, R | X ) exp { g ( t − E − ℓ ; θ ) (cid:21) . Define also Z b ( t ) = A g θ ( t − E − ℓ ) , Z u ( t ) = A g θ ( t − E − ℓ ) + (1 − A ) g θ ( t − R − ℓ ) . Then it is straightforward that the observed-data estimating functions are E Λ b { O ; Λ b ( t ) , θ } = d e N b ( t ) − d Λ b ( t ) e Y b ( t ) , E Λ u { O ; Λ u ( t ) , θ } = d e N u ( t ) − d Λ u ( t ) e Y u ( t ) , E θ { O ; Λ b ( · )Λ u ( · ) , θ } = Z T C Z b ( t ) { d e N b ( t ) − d Λ b ( t ) e Y b ( t ) } + Z L T P Z u ( t ) { d e N u ( t ) − d Λ u ( t ) e Y u ( t ) } . Letting e N bi ( t ), e N ui ( t ), e Y bi ( t ), e Y ui ( t ), Z bi ( t ), and Z ui ( t ) denote evaluation at O i in (1), the foregoingdevelopments lead to the set of observed-data estimating equations n X i =1 { d e N bi ( t ) − d Λ b ( t ) e Y bi ( t ) = 0 , n X i =1 { d e N ui ( t ) − d Λ u ( t ) e Y ui ( t ) = 0 , (27) n X i =1 (cid:20)Z T C Z bi ( t ) { d e N bi ( t ) − d Λ b ( t ) e Y bi ( t ) } + Z L T P Z ui ( t ) { d e N ui ( t ) − d Λ u ( t ) e Y ui ( t ) } (cid:21) = 0 . (28)For fixed θ , the estimators for d Λ b ( t ) and d Λ u ( t ) are the solutions to the equations in (27) given by d b Λ b ( t ) = ( n X i =1 e Y bi ( t ) ) − n X i =1 d e N bi ( t ) , d b Λ u ( t ) = ( n X i =1 e Y ui ( t ) ) − n X i =1 d e N ui ( t ) . (29)Substituting these expressions in (28) yields, after some algebra, the equation n X i =1 (cid:20)Z T C { Z bi ( t ) − Z b ( t ) } d e N bi ( t ) + Z L T P { Z ui ( t ) − Z u ( t ) } d e N ui ( t ) (cid:21) = 0 , (30) Z b ( t ) = ( n X i =1 e Y bi ( t ) ) − n X i =1 Z bi ( t ) e Y bi ( t ) , Z u ( t ) = ( n X i =1 e Y ui ( t ) ) − n X i =1 Z ui ( t ) e Y ui ( t ) . Practical Implementation and Inference
Choice of the weight functions e w ( t, e ), e w ( t, e ), w ( t, e, r ), and w ( t, e, r ) is arbitrary but can play animportant role in performance of the resulting estimators. We recommend taking a fixed value e x of X ,e.g., the sample mean, and setting e w a ( t, e ) = h a ( t, e | e x ) and w a ( t, e, r ) = h a ( e, r | e x ), a = 0 ,
1, wherethe latter does not depend on t . The resulting weights h a ( t, e | e x ) /h j ( t, e | X ) and h a ( e, r | e x ) /h a ( e, r | X ), a = 0 ,
1, are referred to as stabilized weights (Robins, Hern´an, and Brumback, 2000), as they mitigatethe effect of small inverse probability weights that can give undue influence to a few observations. Notethat dependence of the inverse probability weights on p A cancels in construction of stabilized weights.Moreover, if there is no confounding, in that λ R,j ( r | X, A ), j = 1 , f E | X ( e | X ), and p Ψ ( X, Γ) donot depend on X , the stabilized weights are identically equal to one.If the “survival probabilities” for R , K R,j ( r | X, A ), and the densities f R,j ( r | X, A ), j = 1 ,
2, and f E | X ( e | X ) in the inverse probability weights, which appear in the expressions in the estimating equation(30), were known, (30) could be solved to yield an estimator for θ and in particular θ characterizing VEwaning. As these quantities are unknown, models must be posited for them, leading to estimators thatcan be substituted in (30). We propose the use of Cox proportional hazards models for λ R,j ( r | X, A ), j = 1 ,
2, in (19), which can be fitted using the data { X i , A i , R i , I(Γ i = j ) } , i = 1 . . . , n ; and for thehazard of entry time E given X , which can be fitted using ( E i , X i ), i = 1 , . . . , n . A binary, e.g., logistic,regression model can be used to represent p Ψ ( X, Γ) and fitted using ( X i , Γ i , Ψ i ) for i such that A i = 0.For individual i , the stabilized weights involve the quantities f R,j ( R i | e x, a ) /f R,j ( R i | X i , a ), j = 1 , a = 0 ,
1, and f E | X ( E i | e x ) /f E | X ( E i | X i ). With proportional hazards models as above with predictors φ j ( X, β j ), say, it is straightforward that f E | X ( E i | e x ) /f E | X ( E i | X i ) and f R,j ( R i | e x, a ) /f R,j ( R i | X i , a ) = [exp { φ j ( e x, β j ) }K R ( R i | e x, a )] / [exp { φ j ( X i , β j ) }K R ( R i | X i , a )] , where in each case the baseline hazard cancels from numerator and denominator. for Thus, the estimatedstabilized weights involve only the estimated cumulative hazard functions and estimators for the β j , eachof which is root- n consistent and asymptotically normal.As sketched in Appendix F, with stabilized weights set equal to one or estimated, the estimatingequation (30) can be solved easily via a Newton-Raphson algorithm. A heuristic argument demonstratingthat b θ is asymptotically normal leading to an expression for its approximate sampling variance using thesandwich technique is given in Appendix F. 15 Simulations
We report on simulation studies demonstrating performance of the methods, each involving 1000 MonteCarlo replications, based roughly on the Moderna trial. We took p A = 0 . T A = 12, T P = 19, T U = 21, and T C = 31, where all times are in weeks, and consider an analysis at calendar time L = 52weeks, with n = 30,000. In all cases, g ( u, θ ) = θ I( u > v ) where v = 20 weeks and θ = log(0 . v , so that, depending on θ , VE potentially wanes following v . We consider θ = log(7), corresponding to VE = 65% after time v , and θ = 0, corresponding to nowaning.Because the trial and unblinding process are ongoing, we were not able to base our generative scenarioson data from the trial. Owing to the complexity of the trial and multiple potential sources of confounding,to facilitate exploration of a range of conditions while controlling computational complexity and intensity,we focused on several basic scenarios meant to represent varying degrees of confounding consistent withour expectations for the most likely sources of such confounding in the trial. Specifically, we took f E | X ( e | X ) and λ R, ( r | X, A ) to not depend on X (or A in the latter case) in any scenario, reflecting mostlyrandom entry and PDCV unblinding processes. In scenarios involving confounding, we took λ R, ( r | X, A ),corresponding to the period [ T P , T U ) in which “requested unblinding” occurred, and the “agreementprocess” p Ψ ( X, Γ) to depend on X , as described below, reflecting our belief that these processes could beassociated with participant characteristics.In the first set of simulations, we consider two cases: (i) no confounding, where all of λ R,j ( r | X, A ), j = 1 , f E | X ( e | X ), and p Ψ ( X, Γ) do not depend on X ; and (ii) confounding, where λ R, ( r | X, A ) and p Ψ ( X, Γ) depend on X as above. In both (i) and (ii), the entry process E ∼ U (0 , T A ), i.e., uniformon [0 , T A ], and the unblinding process during PDCVs was U ( T U , T C ); see below. In each simulationexperiment, for each participant in each Monte Carlo data set, we first generated A ∼ Bernoulli( p A ), twobaseline covariates X ∼ Bernoulli( p X = 0 .
5) and X ∼ N ( µ X = 45 , σ X = 10 ), and E as above. Toobtain R , we generated G to be exponential with hazard λ R, ( r | X, A ) = exp[ e β + { e β ( X − p X ) + e β ( X − µ X ) } (1 − A ) + { e β ( X − p X ) + e β ( X − µ X ) } A ], where e β = log(0 . T P , T U ), and ( e β , e β , e β , e β ) = (0 , , , ,
0) for (i), no confounding,and ( − . , − . , . , .
08) for (ii), confounding. With R = T P + G and R ∼ U ( T U , T C ), we let e Γ = 1 + I( R ≥ T U ) and e R = R I( e Γ = 1) + R I( e Γ = 2). We generated Ψ as Bernoulli { p Ψ ( X, e Γ) } , p Ψ ( X, e Γ) = expit { e γ + e γ ( X − p X ) + e γ ( X − µ X ) + e γ e Γ } , expit( u ) = (1 + e − u ) − , where e γ = 1 . e γ , e γ , e γ ) = (0 , , − .
1) for (i) and = ( − . , − . , − .
1) for (ii).16o generate U, ∆, we first generated T ∗ ( E, R ) and T ∗ ( E, R ) based on (10)-(11), with λ b ( t ) = λ b =exp { δ + δ ( X − p X ) + δ ( X − µ X ) + Z} , where ( δ , δ , δ ) = { log(0 . , . , . } , leading toapproximately a 3% infection rate for placebo participants over L , and Z ∼ N (0 , . λ uℓ ( t ) = λ uℓ = λ b ; ζ ( t ) = 0; and λ u ( t ) = λ u = 1 . λ b , so that λ a ( t, e, r ) in (10)-(11), a = 0 ,
1, are piecewise constant hazards. T ∗ ( E, R ) and T ∗ ( E, R ) were obtained via inverse transform sampling. We then generated U (calendartime) as U = E + AT ∗ ( E, R ) + (1 − A ) (cid:2) I { T ∗ ( E, R ) < e R } T ∗ ( E, R ) + I { T ∗ ( E, R ) ≥ e R }{ Ψ T ∗ ( E, R ) + (1 − Ψ) T ∗ r } , where T ∗ r = e R + G for G exponential with hazard λ b ; infection times for unblinded placeboparticipants who decline vaccine are not used in the analysis. Finally, we set ∆ = I( U < L ), and defined R = U I( U ≤ e R ) + e R I( U > e R ) and Γ = e ΓI(
U > R ). Although we obtained Ψ for all n participants, Ψ isused only when A = 0, Γ ≥ θ = log(7) and (b) θ = 0, we estimated θ and thus V E ( τ ) for τ ≤ v and τ > v two ways: taking the stabilized weights equal to one, so disregarding possibleconfounding, and with estimated stabilized weights. The latter were obtained by fitting proportionalhazards models for entry time E with linear predictor ν X + ν X and for λ R,j ( r | X, A ), j = 1 , , withlinear predictors β X + β X + β A + β X A + β X A and β X + β X , respectively; and a logisticregression model for p Ψ ( X, Γ) = expit { ( γ + γ X + γ X )I(Γ = 1) + ( γ + γ X + γ X )I(Γ = 2) } .Table 2 presents the results for estimation of θ , dictating waning; V E ≤ = 1 − exp( θ ), VE priorto v = 20 weeks; and V E > = 1 − exp( θ + θ ), VE after v = 20 weeks. Because the Monte Carlodistribution of some of these quantities exhibited slight skewness, those for the VE quantities likely dueto the exponentiation, we report both Monte Carlo mean and median. Estimation of V E ≤ showsvirtually no bias for both (a) and (b); that for V E > in case (a) shows minimal bias and virtually nonefor (b). In all cases, standard errors obtained via the sandwich technique as outlined in Appendix Falong with the delta method for the VEs track the Monte Carlo standard deviations. Under both (i) noconfounding and (ii) confounding, estimation of the stabilized weights appears to have little consequencefor precision of the estimators relative to setting them to equal to one. 95% Wald confidence intervals,exponentiated for the VEs, achieve nominal coverage. For (b) and each combination of stabilized weightsset equal to one or estimated and (i), no confounding, and (ii), confounding, we also calculated theempirical Type I error achieved by a Wald test at level of significance 0.05 for VE waning addressingthe null and alternative hypotheses H : θ ≤ H : θ >
0. These values are 0.043 and 0.056when using stabilized weights set equal to one under (i) and (ii), respectively; the analogous values withestimated weights are 0.046 and 0.050 under (i) and (ii).In the first set of simulations, the confounding induced by our generative choices led to little to no17 able 2: Simulation results based on 1000 Monte Carlo replications, first scenario. Mean = mean of MonteCarlo estimates, Med = median of Monte Carlo estimates, SD = standard deviation of Monte Carlo estimates,SE = average of standard errors obtained via the sandwich technique/delta method, Cov = empirical coverage ofnominal 95% Wald confidence interval (transformed for
V E ). V E ≤ = 1 − exp( θ ) , VE prior to v = 20 weeks; V E > = 1 − exp( θ + θ ) , VE after v = 20 weeks. True values: (a) θ = log(7) = 1 . , V E ≤ = 0 . , V E > = 0 . ; (b) θ = 0 , V E ≤ = V E > = 0 . . Stabilized Weights = 1 Stabilized Weights EstimatedMean Med SD SE Cov Mean Med SD SE Cov(i), no confounding; (a) θ = log(7) θ V E ≤ V E > θ = log(7) θ V E ≤ V E > θ = 0 θ -0.020 -0.019 0.433 0.422 0.954 0.007 0.019 0.421 0.424 0.958 V E ≤ V E > θ = 0 θ V E ≤ V E > θ and the VEs prior to and after 20 weeks. Notably, modeling and fitting of thestabilized weights to adjust for potential confounding shows little effect relative to setting the stabilizedweights to one. To the extent that this scenario is a plausible approximation to actual conditions of thetrial, it may be that confounding will not be a serious challenge for the analysis of VE waning.To examine the ability of the methods with estimated stabilized weights to adjust for confounding thatpotentially could be sufficiently strong to bias results, we carried out additional simulations under settings(a) θ = log(7) and (b) θ = 0 with (ii) confounding in which our choices of generative parameters inducea stronger association between the potential infection times and the agreement process. Specifically, wetook instead ( δ , δ , δ ) T = { log(0 . , . , . } T and ( e γ , e γ , e γ , e γ ) = (1 . , − . , − . , − . able 3: Simulation results based on 1000 Monte Carlo replications, second scenario. Entries are as in Table 2.True values: (a) θ = log(7) = 1 . , V E ≤ = 0 . , V E > = 0 . ; (b) θ = 0 , V E ≤ = V E > = 0 . . Stabilized Weights = 1 Stabilized Weights EstimatedMean Med SD SE Cov Mean Med SD SE Cov(ii), confounding; (a) θ = log(7) θ V E ≤ V E > θ = 0 θ V E ≤ V E > θ and V E > are slightly biased when stabilizedweights are set equal to one, although coverage probability for the latter is at the nominal level. Thisfeature is mitigated by use of estimated stabilized weights. Coverage probability for θ is somewhat lowerthan nominal. Under (b), empirical Type I error achieved by a Wald test at level of significance 0.05 of H : θ ≤ H : θ >
0. is 0.122 when stabilized weights are equal to one, demonstrating thepotential for biased inference; Type I error is 0.065 using estimated stabilized weights, leading to a morereliable test.
We have proposed a conceptual framework based on potential outcomes for study of VE in which as-sumptions on biological, behavioral, and other phenomena are made transparent. The correspondingstatistical framework combines information from blinded and unblinded participants over time. We focuson the setting of the Moderna phase 3 trial, but the principles can be adapted to other settings, includingthe blinded crossover design of Follmann et al. (2020). The methods provide a mechanism to account forpossible confounding.Through condition (ii) in Section 3, (ii) E { π ( t, τ ) | X } /E { π ( t ) | X } = q ( τ ), the methods embed theassumption that VE is similar across current and emerging viral variants. If the analyst is unwilling toadopt an assumption like condition (ii), then it is not possible to rule out that the data from the blinded(prior to T P ) and unblinded (starting at T P phases of the trial reflect very different variant mixtures. Inthis case, calendar time and time since vaccination cannot be disentangled, and thus it is not possible to19valuate VE solely as a function of time since vaccination. However, it may be possible to evaluate theratio of infection rates under vaccine at any time t (and thus variant mixture in force at t ) after differenttimes since vaccination τ and τ , say, during the unblinded phase of the trial, namely, I u ( t, τ ) / I u ( t, τ ), t ≥ T P . The infection rates can be estimated based on the infection status data at time t from vaccinatedindividuals who received vaccine at times t − τ and t − τ , respectively. These infection rates and theirratio will reflect information about the waning of the vaccine itself under the conditions at time t , andin fact this infection rate ratio can be viewed as the ratio of vaccine efficacies at τ and τ . However,because after T C information on I u ( t ) will no longer be available, it is not possible to deduce VE itselffor t ≥ T C . But if data external to the trial became available that provide information on VE at t , evenfor small τ , it may be possible to integrate this information with that from the infection rates to gaininsight into VE as a function of τ . Acknowledgements
The authors thank Dean Follmann for helpful discussions and insights.
References
Baden, L. R., El Sahly, H. M, Essink, B., Kotloff, K., Frey, S., Novak, R., et al. for the COVE StudyGroup (2020). Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine.
New England Journalof Medicine , https://doi.org/10.1056/NEJMoa2035389 .Fintzi, J. and Follmann, D. (2021). Assessing vaccine durability in randomized trials following placebocrossover. arXiv preprint arXiv:2101.01295v2.Fleming, T. R. and Harrington, D. P. (2005). Counting Processes and Survival Analysis . New York:Wiley.Follmann, D., Fintzi, J., Fay, M. P., Janes, H. E., Baden, L., El Sahly, H. et al. (2020). Assessingdurability of vaccine effect following blinded crossover in COVID-19 vaccine efficacy trials. medRxiv.2020 Dec 14;2020.12.14.20248137. https://doi.org/10.1101/2020.12.14.20248137 .Halloran, M. E., Longini, I. M, and Struchiner, C. J. (1996). Estimability and interpretation of vaccineefficacy using frailty mixing models.
American Journal of Epidemiology , 144, 83–97.20in, D.-Y., Zeng, D., and Gilbert, P. B. (2021). Evaluating the long-term efficacy of COVID-19 vaccines.medRxiv. 2021 Jan 13;2021.01.13.21249779. https://doi.org/10.1101/2021.01.13.21249779 .Longini, I. M. and Halloran, M. E. (1996). frailty mixture model for estimating vaccine efficacy.
Journalof the Royal Statistical Society, Series C , 45, 165–173.Moderna Clinical Study Protocol, Amendment 6, 23 December 2020, available at
Robins, J. M., Hern´an, M. A., and Brumback, B. (2000). Marginal structural models and causal inferencein epidemiology.
Epidemiology , 11, 550–560.Rubin, D. B. (1980). Bias reduction using Mahalanobis-metric matching.
Biometrics , 36, 293–298.Yang, S., Tsiatis, A. A., and Blazing, M. (2018). Modeling survival distribution as a function of timeto treatment discontinuation: A dynamic treatment regime approach.
Biometrics , 74, 900–909.
Appendix A: Demonstration of (2) and (3)
We demonstrate that, under the conditions in Section 3 of the main paper, namely,(i) { π ( t, τ ) , π ( t ) } ⊥⊥ { S, c b ( t ) }| X and { π ( t, τ ) , π ( t ) } ⊥⊥ { S, c u ℓ ( t ) , c u ( t ) }| X ,(ii) E { π ( t, τ ) | X } /E { π ( t ) | X } = q ( τ ),that (2) of the main paper, R b ( t, τ ) = I b ( t, τ ) I b ( t ) = E { p ( t, S ) c b ( t ) π ( t, τ ) } E { p ( t, S ) c b ( t ) π ( t ) } (A.1)does not depend on t , and the second equality in (3) of the main paper, I u ( t, τ ) I u ( t, τ ) = E { p ( t, S ) c u ( t ) π ( t, τ ) } E { p ( t, S ) c u ( t ) π ( t, τ ) } = R b ( τ ) R b ( τ ) , τ , τ ≥ ℓ ; (A.2)the first equality in (3) of the main paper follows by an entirely similar argument.We can write (A.1) using condition (i) as R b ( t, τ ) = E (cid:2) E { p ( t, S ) c b ( t ) | X } E { π ( t, τ ) | X } (cid:3) E (cid:2) E { p ( t, S ) c b ( t ) | X } E { π ( t ) | X } (cid:3) .
21y condition (ii), E { π ( t, τ ) | X } = E { π ( t ) | X } q ( τ ); thus, substituting yields R b ( t, τ ) = E (cid:2) E { p ( t, S ) c b ( t ) | X } E { π ( t ) | X } (cid:3) q ( τ ) E (cid:2) E { p ( t, S ) c b ( t ) | X } E { π ( t ) | X } (cid:3) = q ( τ ) , so that in fact q ( τ ) = R b ( τ ).We can write (A.2) as E (cid:2) E { p ( t, S ) c u ( t ) | X } E { π ( t, τ ) | X } (cid:3) E (cid:2) E { p ( t, S ) c u ( t ) | X } E { π ( t, τ ) | X } (cid:3) . Under condition (ii), E { π ( t, τ j ) | X } = E { π ( t ) | X } q ( τ j ), j = 1 ,
2; thus, substituting these equalitiesyields E (cid:2) E { p ( t, S ) c u ( t ) | X } E { π ( t ) | X } q ( τ ) (cid:3) E (cid:2) E { p ( t, S ) c u ( t ) | X } E { π ( t ) | X } q ( τ ) } (cid:3) = q ( τ ) q ( τ ) = R b ( τ ) R b ( τ ) , as required. Appendix B: Discussion of Assumptions
The conceptual framework in Section 3 of the main paper in which we define vaccine efficacy at a particulartime since vaccination relies on some assumptions. Of critical importance is the assumption referred toas condition (ii), namely, E { π ( t, τ ) | X } /E { π ( t ) | X } = q ( τ ) , (B.1)which states that the ratio of transmission probabilities over time within values of X does not changewith time and does not depend on characteristics in X but depends only on time since vaccination.As noted in Section 3 of the main paper, in our conceptualization, we let the individual-specifictransmission probabilities π ( t, τ ) and π ( t ) depend on t to reflect an evolving mixture of viral variantsas mutations of the virus occur over the course of the pandemic, under which the overall virulence ofvirus to which individuals in the study population may be exposed is changing. From this point ofview, we can regard time t as a “proxy” for this changing variant mixture and its virulence as the studyprogresses. If in fact the overall virulence of the variant mixture does not change or changes only graduallyover time, then it may be reasonable to take π ( t, τ ) = π ( τ ) and π ( t ) = π . In this case, the ratio E { π ( t, τ ) | X } /E { π ( t ) | X } in (B.1) is a function only of τ and X . If instead the variant mixture doeschange over the course of the study in a non-trivial way, taking π ( t, τ ) and π ( t ) not to depend on t isuntenable. However, if within the mixture of variants present at any time t we are willing to assume thatthe ratio of transmission probabilities between vaccine and placebo stays in constant proportion for allvariants, it again is reasonable to assume that E { π ( t, τ ) | X } /E { π ( t ) | X } does not depend on t so is afunction only of τ and X . 22nder either of these perspectives, for (B.1) to hold, we furthermore must be willing to assume that E { π ( t, τ ) | X } /E { π ( t ) | X } does not depend on X (in addition to not depending on t ) and thus dependsonly on τ . Adopting (B.1) is similar in spirit to making the assumptions embodied in many popularmodels; e.g., a constant odds ratio over categories in the proportional odds model or a constant hazardratio over time in the proportional hazards model. If (B.1) is violated in that E { π ( t, τ ) | X } /E { π ( t ) | X } does depend on X (but not on t ), the implication for the proposed methods is that, in estimating VEassuming it depends only on τ , one is estimating roughly a weighted average of VE as a function of τ overvalues of X in a manner similar to the Mantel-Haenzel method; such an interpretation is also commonlyinvoked when the proportional odds or hazards assumptions do not hold.If the analyst is unwilling to adopt an assumption like that in (B.1), then it is not possible to rule outthat the data from the blinded (prior to T P ) and unblinded (starting at T P , when unblinding requestscommenced following the Pfizer EUA) phases of the trial reflect very different variant mixtures. In thiscase, calendar time and time since vaccination cannot be disentangled, and thus it is not possible toevaluate vaccine efficacy solely as a function of time since vaccination. In this setting, however, it maybe possible to evaluate the ratio of infection rates under vaccine at any time t (and thus variant mixturein force at t ) after different times since vaccination τ ≥ ℓ and τ ≥ ℓ , say, during the unblinded phase ofthe trial, namely, I u ( t, τ ) / I u ( t, τ ) , t ≥ T P . The infection rates in this ratio presumably can be estimated based on the infection status data at time t from vaccinated individuals who received vaccine at times t − τ and t − τ , respectively. These infectionrates and their ratio will reflect information about the waning of the vaccine itself under the conditionsat time t , and in fact this infection rate ratio can be viewed as the ratio of vaccine efficacies at differentvalues τ and τ . However, because after T C information on I u ( t ) will no longer be available, it is notpossible to deduce vaccine efficacy itself for t ≥ T C . But if data external to the trial became availablethat provide information on vaccine efficacy at t , even for small τ it may be possible to integrate thisinformation with that from the infection rates to gain insight into vaccine efficacy itself as a function of τ . 23 ppendix C: Approximate Equivalence of Hazard Rate and InfectionRate As an example, consider λ ( t, e ) = λ ( t, e, ∞ ) defined in (9) of the main paper. From Section 3 of themain paper, the individual-specific infection rate for an arbitrary subject in the study population at site S who receives placebo and is never unblinded ( r = ∞ ) is given by p ( t, S ) c b ( t ) π ( t ). This quantity isa random variable defined for the population Ω with probability { P ( ω ) : ω ∈ Ω } , where we view ω asan individual in Ω. Thus, the infection rate for ω ∈ Ω is ι ( t )( ω ) = p { t, S ( ω ) } c b ( t )( ω ) π ( t )( ω ), and thepopulation-level infection rate is given by E { ι ( t ) } = Z Ω ι ( t )( ω ) dP ( ω ) . In contrast, the hazard at time t is defined by λ ( t, e ) = − ddt log (cid:2) pr { T ∗ ( e ) + e ≥ t } (cid:3) , where pr { T ∗ ( e ) + e ≥ t } = Z Ω pr { T ∗ ( e )( ω ) + e ≥ t } dP ( ω ) . If ω is at risk of infection at time t , then this individual’s hazard of becoming infected at t is given by ι ( t )( ω ). Thus, pr { T ∗ ( e )( ω ) + e ≥ t } = exp (cid:26) − Z te ι ( u )( ω ) du (cid:27) . We make the rare infection assumption Z L ι ( t )( ω ) du < ǫ a.s. (C.1)Now λ ( t, e ) = R Ω G ( t )( ω ) ι ( t )( ω ) dP ( ω ) R Ω G ( t )( ω ) dP ( ω ) , where, using the rare infection assumption, G ( t )( ω ) = exp (cid:26) − Z te ι ( u )( ω ) du (cid:27) ≥ exp( − ǫ ) > − ǫ a.s.Because Z Ω G ( t )( ω ) ι ( t )( ω ) dP ( ω ) ≤ Z Ω ι ( t )( ω ) dP ( ω )and G ( t )( ω ) > − ǫ a.s., λ ( t, e ) ≤ R Ω ι ( t )( ω ) dP ( ω )1 − ǫ . Z Ω G ( t )( ω ) ι ( t )( ω ) dP ( ω ) > (1 − ǫ ) Z Ω ι ( t )( ω ) dP ( ω )and R Ω G ( t )( ω ) dP ( ω ) ≤ λ ( t, e ) ≥ (1 − ǫ ) Z Ω ι ( t )( ω ) dP ( ω ) . Thus, (1 − ǫ ) < λ ( t, e ) R Ω ι ( t )( ω ) dP ( ω ) < (1 − ǫ ) − . Consequently, under the rare infection assumption (C.1), the population-level infection rate and thepopulation-level hazard rate are of the same order of magnitude.
Appendix D: Derivation of Estimating Functions (12)-(14)
We present derivations leading to the estimating functions (12)-(14) based on potential outcomes given inSection 4.2 of the main paper. Because interest focuses on τ ≥ ℓ , from (10) and (11) of the main paper, weare concerned only with Λ b ( t ), Λ u ( t ), and θ . Accordingly, to determine appropriate linear combinationsof the mean-zero counting process increments { dN ∗ a ( t, e, r ) − d Λ a ( t, e, r ) Y ∗ a ( t, e, r ) } , a = 0 ,
1, we mustdeduce relevant values of t , e , and r , where e ≤ T A and T P ≤ r < T C by design. For a = 0, from (10) ofthe main paper, the relevant values are t < r or t ≥ ℓ + r and e ≤ min( t, r ). For a = 1, from (11) of themain paper, e + ℓ ≤ t ≤ r and t > r .Consider for fixed 0 ≤ t ≤ L , a = 0 ,
1, integrals of the form
Z Z { dN ∗ a ( t, e, r ) − d Λ a ( t, e, r ) Y ∗ a ( t, e, r ) } w a ( t, e, r ) dr de, (D.1)where w a ( t, e, r ) is a non-negative weight function, a = 0 ,
1. We determine the limits of integration for(D.1) by considering three time periods.When t < T P , at which point all trial participants are still blinded, so that t < r , (D.1) for a = 0becomes, using (10) of the main paper and the consistency assumptions below (11) of the main paper, Z min( t, T A )0 Z T C T P { dN ∗ ( t, e ) − d Λ b ( t ) Y ∗ ( t, e ) } w ( t, e, r ) dr de = Z min( t, T A )0 { dN ∗ ( t, e ) − d Λ b ( t ) Y ∗ ( t, e ) } e w ( t, e ) de, (D.2)where for t < T P e w ( t, e ) = Z T C T P w ( t, e, r ) dr. a = 1, ℓ ≤ t < T P shows that (D.1) becomes, using (11) of the main paper, Z min( t − ℓ, T A )0 Z T C T P (cid:2) dN ∗ ( t, e ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ ) } Y ∗ ( t, e ) (cid:3) w ( t, e, r ) dr de = Z min( t − ℓ, T A )0 (cid:2) dN ∗ ( t, e ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ ) } Y ∗ ( t, e ) (cid:3) e w ( t, e ) de, (D.3)where for t < T P e w ( t, e ) = Z T C T P w ( t, e, r ) dr. Next consider T P ≤ t < T C ; at times in this interval, some participants are still blinded while othershave become unblinded. We consider both t < r , so before unblinding, and t ≥ r , after unblinding attime r . First consider (D.1) with a = 0. For t < r , (D.1) becomes Z T A Z T C t { dN ∗ ( t, e ) − d Λ b ( t ) Y ∗ ( t, e ) } w ( t, e, r ) dr de = Z T A { dN ∗ ( t, e ) − d Λ b ( t ) Y ∗ ( t, e ) } e w ( t, e ) de, (D.4)where for T P ≤ t < T C e w ( t, e ) = Z T C t w ( t, e, r ) dr. Similarly, for a = 1, t < r , (D.1) becomes Z T A (cid:2) dN ∗ ( t, e ) − d Λ b ( t ) exp { θ + g ( t − e − ℓ ; θ ) } Y ∗ ( t, e ) (cid:3) e w ( t, e ) de, (D.5)where for T P ≤ t < T C e w ( t, e ) = Z T C t w ( t, e, r ) dr. Continuing to consider T P + ℓ ≤ t < T C , now take t ≥ r . For a = 0, (D.1) becomes Z T A Z t − ℓ T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − r − ℓ ; θ )I( t − r ≥ ℓ ) } Y ∗ ( t, e, r ) (cid:3) w ( t, e, r ) dr de, (D.6)For a = 1, (D.1) becomes Z T A Z t T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − e − ℓ ; θ ) } Y ∗ ( t, e, r ) (cid:3) w ( t, e, r )I( t ≥ r ) dr de, (D.7)Finally, consider t ≥ T C ; these are times where all participants are unblinded. Thus, when a = 0,(D.1) equals Z T A Z min(( t − ℓ, T C ) T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − r − ℓ ; θ )I( t − r ≥ ℓ ) } Y ∗ ( t, e, r ) (cid:3) w ( t, e, r ) dr de, (D.8)26nd when a = 1 equals Z T A Z T C T P (cid:2) dN ∗ ( t, e, r ) − d Λ u ( t ) exp { g ( t − e − ℓ ; θ ) } Y ∗ ( t, e, r ) (cid:3) w ( t, e, r )I( t ≥ r ) dr de. (D.9)Combining (D.2)-(D.5) yields estimating function E Λ b { W ∗ ; Λ b ( t ) , θ } in (12) of the main paper.Combining (D.6)-(D.9) yields E Λ u { W ∗ ; Λ u ( t ) , θ } in (13) of the main paper. Estimating function E θ { W ∗ ; Λ b ( · )Λ u ( · ) , θ } arises through similar considerations, integrating over t and differentiating withrespect to θ and θ . Appendix E: Demonstration of (21)-(24)
We make the assumptions (16)-(20) in Section 4.3 of the main paper. Here, we show the first equalitiesin (21) and (23), i.e., E (cid:26) I ( t, e ) dN ( t ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e ) (E.1)and E (cid:26) I ( t, e, r ) dN ( t ) h ( e, r | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e, r ) . (E.2)Demonstration of the other equalities in (21)-(24) follows by analogous arguments.We first show (E.1). By the consistency assumption (16) in the main paper, the left hand side of(E.1) is equal to E (cid:26) I ( t, e ) dN ∗ ( t, e ) h ( t, e | X ) (cid:12)(cid:12)(cid:12)(cid:12) X, W ∗ (cid:27) = dN ∗ ( t, e ) h ( t, e | X ) E { I ( t, e ) | X, dN ∗ ( t, e ) = 1 , W ∗ } . The result follows if we show that E { I ( t, e ) | X, dN ∗ ( t, e ) = 1 , W ∗ } = h ( t, e | X ) . (E.3)By (15) of the main paper, the left hand side of (E.3) is computed aspr { E = e | X, dN ∗ ( t, e ) = 1 , W ∗ } (E.4) × pr { A = 0 | X, E, dN ∗ ( t, e ) = 1 , W ∗ } (E.5) × pr { R > t | X, A = 0 , dN ∗ ( t, e ) = 1 , W ∗ } , (E.6)where we have used the assumption discussed above (19) in the main paper in (E.6). By (17) of themain paper, (E.5) is equal to pr( A = 0) = 1 − p A . By (17) and (18) of the main paper, (E.4) is equal to f E | X ( e | X ). The proof will be complete by showing that (E.6) is equal to K R ( t | X, A = 0).To demonstrate this, we consider t < T P , T p ≤ t < T U , T U ≤ t < T C , and t ≥ T C in turn. Clearly(E.6) is equal to 1 for t < T P . Because the estimating function using I ( t, e ) is defined only for t < T C ,27e need not consider t ≥ T C . Thus, we need only consider the cases T p ≤ t < T U and T U ≤ t < T C . For T p ≤ t < T U , we write (E.6) as a product integral as in Anderson et al. (1993) and use an argumentsimilar to that in (8.72)-(8.77) of Tsiatis et al. (2020): Y T P ≤ w 4) converges in probability tozero. Thus, solving (F.2) is asymptotically equivalent to setting (F.3) equal to zero. Letting θ (0) denotethe true value of θ under the assumption that the semiparametric model (6) of the main paper is correctlyspecified, then (F.3) is a sum of mean-zero independent and identically distributed (iid) terms ψ ( O i ; θ ) = Z T C { Z bi ( t ) − µ b ( t ) }{ d e N bi ( t ) − d Λ b ( t ) e Y bi ( t ) } + Z L T P { Z ui ( t ) − µ u ( t ) }{ d e N ui ( t ) − d Λ u ( t ) e Y ui ( t ) } , (F.5)30here E { ψ ( O i ; θ (0) ) } = 0. Thus, the estimator b θ solving the asymptotically equivalent estimating equa-tion n X i =1 ψ ( O i ; θ ) = 0satisfies, by a standard Taylor series expansion,0 = n X i =1 ψ ( O i , b θ ) ≈ n X i =1 ψ ( O i , θ (0) ) + ( n X i =1 ∂ψ ( O i , θ ) ∂θ T ) ( b θ − θ (0) ) . As a consequence, n / ( b θ − θ (0) ) = " − E ( ∂ψ ( O i , θ (0) ) ∂θ T ) − n − / n X i =1 ψ ( O i , θ (0) ) + o P (1) , (F.6)which implies that b θ is asymptotically normal with mean zero and covariance matrix " − E ( ∂ψ ( O i , θ (0) ) ∂θ T ) − var { ψ ( O i , θ (0) ) } " − E ( ∂ψ ( O i , θ (0) ) ∂θ T ) − T , (F.7)where var { ψ ( O i , θ (0) ) } = E { ψ ( O i , θ (0) ) ψ ( O i , θ (0) ) T } .An estimator for the asymptotic variance (F.7) can be obtained as follows. The term var { ψ ( O i , θ (0) ) } can be estimated by c var { ψ ( O i , θ (0) ) } = n − n X i =1 b ψ i ( b θ ) b ψ i ( b θ ) T , where b ψ i ( b θ ) is an estimator for ψ ( O i , θ (0) ) obtained by substituting (i) Z k ( t ) for µ k ( t ), k = b, u ; (ii) d b Λ k ( t )in (29) of the main paper for d Λ k ( t ), k = b, u ; and (iii) b θ for θ (0) . An estimator for E ( ∂ψ ( O i , θ (0) ) ∂θ T ) is obtained by substitutions (i)–(iii) in this expression and averaging over i , leading to b E ( ∂ψ ( O i , θ (0) ) ∂θ T ) = − n − n X i =1 (cid:26)Z T C V b ( t ) d e N bi ( t ) + Z L T P V u ( t ) d e N ui ( t ) (cid:27) , (F.8)where V k ( t ) = P ni =1 { Z ki ( t ) − Z k ( t ) }{ Z ki ( t ) − Z k ( t ) } T e Y ki ( t ) P ni =1 e Y ki ( t ) , k = b, u. The resulting sandwich estimator for the large sample covariance matrix of b θ is then given by " b E ( ∂ψ ( O i , θ (0) ) ∂θ T ) − c var { ψ ( O i , θ (0) ) } " b E ( ∂ψ ( O i , θ (0) ) ∂θ T ) − . (F.9)31he foregoing developments take the inverse probability weights and thus the stabilized weights tobe known. If models for λ R,j ( r | X, A ), j = 1 , f E | X ( ( e | X ), and p Ψ ( X, Γ) are posited and fitted andsubstituted in (F.1), then the large sample distribution of n / ( b θ − θ (0) ) would be considerably morecomplicated. In simulations, we have observed that standard errors and confidence intervals based on(F.6) and (F.9) reflect the true sampling variation in that their numerical values are consistent withthe Monte Carlo sampling variation and confidence intervals achieve the nominal level of coverage. Analternative strategy to obtaining approximate standard errors and confidence intervals would be to use anonparametric bootstrap.The result (F.6) suggests a Newton-Raphson iterative scheme for solving the estimating equation(F.1). Letting θ (0) be an initial value for θ and θ ( m ) be the value at the m th iteration, compute theupdate by θ ( m +1) = θ ( m ) − (cid:20) b E (cid:26) ∂ψ ( O i , θ ( m ) ) ∂θ T (cid:27)(cid:21) − b ψ i ( θ ( m ) ) . This scheme is iterated until some convergence criterion is satisfied. Appendix References Anderson, P. K., Borgan, Ø., Gill, R. D., and Keiding, N. (2993). Statistical Methods Based on CountingProcesses . New York: Springer.Tsiatis, A. A., Davidian, M., Holloway, S. T., and Laber, E. B. (2020).