Assessing the prior event rate ratio method via probabilistic bias analysis on a Bayesian network

Background: Unmeasured confounders are commonplace in observational studies conducted using real‐world data. Prior event rate ratio (PERR) adjustment is a technique shown to perform well in addressing such confounding. However, it has been demonstrated that, in some circumstances, the PERR method actually increases rather than decreases bias. In this work, we seek to better understand the robustness of PERR adjustment. Methods: We begin with a Bayesian network representation of a generalized observational study, which is subject to unmeasured confounding. Previous work evaluating PERR performance used Monte Carlo simulation to calculate joint probabilities of interest within the study population. Here, we instead use a Bayesian networks framework. Results: Using this streamlined analytic approach, we are able to conduct probabilistic bias analysis (PBA) using large numbers of combinations of parameters and thus obtain a comprehensive picture of PERR performance. We apply our methodology to a recent study that used the PERR in evaluating elderly‐specific high‐dose (HD) influenza vaccine in the US Veterans Affairs population. That study obtained an HD relative effectiveness of 25% (95% CI: 2%‐43%) against influenza‐ and pneumonia‐associated hospitalization, relative to standard‐dose influenza vaccine. In this instance, we find that the PERR‐adjusted result is more like to underestimate rather than to overestimate the relative effectiveness of the intervention. Conclusions: Although the PERR is a powerful tool for mitigating the effects of unmeasured confounders, it is not infallible. Here, we develop some general guidance for when a PERR approach is appropriate and when PBA is a safer option.


INTRODUCTION
From healthcare to econometrics to the social sciences, the findings of observational studies almost always suffer from bias due to confounding. Numerous methods to control for confounding have been developed, and for each of these, dealing with unmeasured confounders usually presents the trickiest challenge.
One method that has been used to address confounding due to unmeasured variables is prior event rate ratio (PERR) adjustment, developed by Weiner and colleagues and shown to perform well in reproducing the results of a Scandinavian randomized controlled trial (RCT) with an observational study conducted using UK Electronic Medical Record (EMR) data. 1 Subsequent work by Weiner et al showed a similarly good performance of PERR applied to observational studies in replicating other RCTs. [2][3][4] Uddin et al 5 studied the performance of PERR and identified situations in which PERR adjustment increases rather than decreases bias. Further development and refinement, including derivation from PERR adjustment of a pairwise Cox likelihood function, was carried out by Lin and Henley. 6 Here, we endeavor to understand when PERR yields meaningful bounds on the size and/or direction of the treatment effect of interest, even if it cannot be guaranteed to remove (or even reduce) bias. We begin with a Bayesian network 7,8 representation of the causal pathways in a standard observational study ( Figure 1). We use Bayesian network analysis to define exact analytic expressions for the joint probabilities of study variables. (Uddin et al used Monte Carlo simulation to approximate the joint probabilities.) We examine the behavior of the system and identify circumstances where the PERR method overestimates and underestimates the true effect of a treatment. As an illustrative example, we then apply our approach to performing a probabilistic bias analysis (PBA) of a study of the relative effectiveness of high-dose (HD) versus standard-dose (SD) seasonal influenza vaccine in the Veterans Affairs (VA) patient population. 9

THE PERR METHOD
Let us suppose we are conducting a cohort study in which the rate of the outcome of interest is measured not just during the observation period after the treatment but also during a baseline period before the treatment. We denote a possible event occurring during the baseline period as Y 1 , a possible event occurring during the observation period as Y 2 and a possible treatment as X. The variables are dichotomous, so, for example, an individual with ( experienced an event in the baseline period, did not receive the treatment and did not experience an event in the observation period. To streamline the notation, we will write Y 1 = y 1 as y 1 , etc. We use P to denote the probability, over a time interval T, that an event occurs in an individual. For example, P(y 2 | x) is the probability per time T that an individual experiences an event in the observation period, conditional on having received treatment. Measured at the population level, P becomes the average incidence rate, ie, events per person per time T. We denote the incidence rate as R. Since the two quantities are numerically identical (assuming that both are defined over the same unit of time, and that individuals are independent), we will use them interchangeably. The incidence rate ratio (RR) in the baseline period is then whereas, in the observation period, it is Now, let us suppose that we find RR base for the event to differ from unity by a statistically significant amount. This suggests that the treatment and control arms are unbalanced with respect to the distribution of one or more determinants of the event. Assuming that there are no systematic measurement errors and that all measurable confounders have already been FIGURE 1 Bayesian network, or probabilistic directed acyclic graph, with dichotomous variables denoting baseline period event Y 1 , treatment X, observation-period event Y 2 , time-varying unobserved confounder {C 1a , C 1b , C 2 }, and the rate ratios (RRs) describing the associations among them. The RRs are defined in Equations (10) to (18) [Colour figure can be viewed at wileyonlinelibrary.com] controlled for, we are led to suspect that there is residual unmeasured confounding by indication. The PERR method attempts to correct for the imbalance and thus recover the "true" RR of the treatment through normalizing RR obs by RR base

THE CAUSAL MODEL
We begin with a Bayesian network (or probabilistic directed acyclic graph [DAG]) 7 depicting, from the perspective of a single individual, the potential causal associations in our study. We include an unmeasured dichotomous confounder, We define distinct values for the confounder in the baseline period (C 1 ) and the observation period (C 2 ). We further divide C 1 into C 1a and C 1b to allow for the possibility that the relationship between C 1 and Y 1 is bidirectional (ie, the state of C 1 influences the state of Y 1 , and the state of Y 1 influences the state of C 1 ). We also allow for the possibility of a direct causal connection between baseline event Y 1 and treatment X. One can show (see the work of Greenland et al 7 ) that this DAG has a null adjustment set and that there exist no instruments or conditional instruments that might permit instrumental variable regression.
We write down a set of equations that describe this causal diagram. The effect of each directed edge is expressed as an RR whose value depends on the state of the variable corresponding to the edge's originating vertex. The model equations describing the population-level incidence rates R (or equivalently, individual-level probabilities P) over a time period T of the occurrence of c, y 1 , x, and y 2 are where r c1a , r y1 , r c1b , r x , r c2 , r Y2 are constants on [0,1]. The operator Π [0,1] (x), ∈ R, is the closest-element mapping from real numbers to [0,1], ie, The RRs are defined as follows: where the F's are constants >0. Uddin et al used Monte Carlo simulation to obtain approximate conditional rates/probabilities from this model. Instead, we will use Bayesian network analysis (see, eg, the work of Pearl 8 ) to calculate exact probabilities, as follows.
The incidence rate with which a given set of values of C 1a , Y 1 , C 1b , X, C 2 , and Y 2 is realized on the above network is their joint probability Joint rates/probabilities of subsets of the variables are calculated by summing the joint probabilities over all possible combinations of the remaining variables, for example, Conditional rates/probabilities can likewise be calculated, for example, We implement these calculations in the R language. 10

THE PERFORMANCE OF PERR ADJUSTMENT
In an observational study, the following incidence rates will be measured: • R(x), the rate of treatment across the study population; • R( 2 |x), the rate of the event in the treatment arm during the observation period; , the rate of the event in the control arm during the observation period; • R( 2 ), the rate of the event across the whole study population during the observation period (where R( 2 ) = R( 2 |x)R(x) + R( 2 |x)R(x)).
If patients have also been observed during a baseline period, then we further have measurements of the following: • R( 1 |x), the rate of the event in the treatment arm during the baseline period (ie, the rate among those who are later treated); • R( 1 |x), the rate of the event in the control arm during the baseline period (ie, the rate among those who are not later treated); • R( 1 ), the rate of the event across the whole study population during the baseline period (where R( 1 ) = R( 1 |x)R(x) + R( 1 |x)R(x)).

FIGURE 2
The conventional rate ratio (RR) (RR obs ), the prior event rate ratio (PERR) estimator, and the true RR (=F x − Y 2 , set to 0.7) in the case where The RRs between treatment and control arms in the observation and baseline period are given by Equations (1) and (2), respectively, while the PERR estimate of the treatment effect is given by Equation (3). The true, unbiased effect of treatment x on the probability of event y 2 is given by F x → Y2 (see Equation (18)). The question we wish to answer for our causal model is, under what circumstances does PERR succeed in decreasing bias, and under what circumstances does it actually increase bias? Furthermore, under what circumstances is the PERR estimate either <F x → Y2 or >F x → Y2 ? For specificity, we take events y to be harmful, thus the smaller F x → Y2 , the more effective the treatment x. Therefore, when PERR/F x → Y2 < 1, the effectiveness of x is overestimated (overoptimistic PERR), whereas if PERR/F x → Y2 > 1, it is underestimated (pessimistic PERR).
We begin by computing the behavior of the PERR in comparison to the true effect F x → Y2 under different scenarios.
1. Unobserved confounder effect is present and affects the baseline and observation period event probability equally (ie, F c2 → Y2 = F c1a → Y1 ≠ 1) and also affects the probability of treatment (ie, F c1b → X ≠ 1), while F y1 → X is fixed at 1. Baseline period events do not affect the probability that the confounder is present (ie, F y1 → C1b = 1). Figure 2 shows RR obs , the PERR, and true effect F x → Y2 as a function of F c1a → Y1 . As we can see, even though RR obs diverges from the true effect, the PERR is exactly equal to the true effect. 2. Retaining the effect of the unobserved confounder, we allow F y1 → Y2 , the effect of the occurrence of a baseline event on the likelihood of an observation period event, to vary. Results are shown in Figure 3. As can be seen, increasing F y1 → Y2 above 1 makes the PERR progressively more pessimistic (ie, PERR/F x → Y2 increases above 1), though the effect is relatively weak and bounded. 3. We allow F c2 → Y2 to differ from F c1a → Y1 . Results are shown in Figure 4. When F c1a → Y1 < F c2 → Y2 , (ie, the effect of the confounder on likelihood of an event is weaker in the baseline period than in the observation period), the PERR is pessimistic; when the converse is true, F c1a → Y1 > F c → Y2 , then the PERR is overoptimistic. As we will see later

FIGURE 3
Performance of the prior event rate ratio (PERR) estimator as a function of F y1 → Y 2 . The ratio of the PERR to the true effect, ie, F x → Y 2 , is plotted. A ratio greater than 1 thus corresponds to an overestimate of the true RR and, hence, an underestimate of the true treatment effect. Scenarios with F c1a → Y 1 = 1 (Scenario a), = 2 (Scenario b), = 3 (Scenario), and = 4 (Scenario d) are shown. As before, we set C 1a = C 1b = C 2 , F c2 → Y 2 = F c1a → Y 1 , and fix F c1b → X = 1.5, r c1a = 0.25. In each case, the overestimate of the true rate ratio (RR) by the PERR initially becomes greater with F y1 → Y 2 , then reaches a maximum, and decreases again. The height of the maximum increases with F c1a → Y 1 . When F c1a → Y 1 = 4, the maximum PERR overestimate of the true RR is 15% and occurs when F y1 → Y 2 = 4

FIGURE 4
Accuracy of the prior event rate ratio (PERR) estimator as a function of F c2 → Y 2 , with F c1a → Y 1 held fixed at 3. In Scenario a, F c1b → X = 1.5, whereas in Scenario b, F c1b → X = 1.9. As before, r c1a = 0.25

FIGURE 5
Performance of the prior event rate ratio (PERR) estimator when F y1 → X differs from 1. In Scenario a, there are no confounders, ie, r c1a = 0, and as before, C 2 = C 1b = C 1a . In Scenario b, the confounder is present (Section 5 and Appendix A), one reason that these two values may differ is if the effect of treatment x in an individual changes depending on whether the confounder is present or absent. 4. We allow F y1 → X to differ from 1; see Figure 5. Above 1, even modest values of F y1 → X suffice to make the PERR substantially overoptimistic. Conversely, as F y1 → X is decreased below 1, PERR rapidly becomes more pessimistic. Note that F y1 → X cannot exceed r −1 X , otherwise P(X = x) would exceed 1 for individuals having a baseline event y 1 . 5. We reverse the directionality of the baseline period relationship between confounder and event: while a pre-existing confounder c 1 does not affect the probability of a baseline event y 1 , an individual without pre-existing c 1 may develop c 1 as a result of experiencing y 1 . Figure 6 shows that, as the strength of the association F y1 → C1b increases, PERR becomes increasingly overoptimistic. Depending on the pre-existing confounder incidence rate r c1a , and on F c1b → X , PERR/F x → Y2 may or may not have a lower limit >0. 6. We allow for the possibility that an individual either loses or gains the confounder between baseline and observation period. We see from Figure 7 that both scenarios act to make PERR overoptimistic; by how much depends on F c1a → Y1 and F c1b → X .
We thus see that the PERR works just as intended as long as the association between confounder and events has the same strength in the baseline and observation periods. The more confounding differs between the periods, though, the poorer the PERR does as an estimator of the true effect. Whether the strength of confounding increases or decreases with time is important, because this determines whether the PERR overestimates or underestimates, respectively, the true effect of the intervention.
We also see that the PERR is overoptimistic when there is a causal connection directed from baseline period event y 1 to treatment x, either directly or via confounder C 1b . Finally, PERR is overoptimistic when individuals are able to either gain or lose the confounder over the course of the study period.

FIGURE 6
Performance of the prior event rate ratio (PERR) estimator when the directionality of the causal relationship between baseline confounder and Y 1 is reversed, ie, F c1a → Y 1 = 1 (presence of c 1a does not affect probability that y 1 is present), F y1 → C1b > 1 (presence of y 1 increases the probability that c 1b is present), and as above, c 1b = c 1a . Scenario a: r c1a = 0.25. Scenario b: r c1a = 0.5. Scenario c: r c1a = 0.25, F c1b → X is increased from 1.5 to 3

Probabilistic bias analysis
We have investigated ways in which using the PERR to control for unobserved confounding can fail. However, if one has sufficient information to be able to place limits on the strengths of the causal associations among the study subjects' set of attributes {C 1a , Y 1 , C 1b , X, C 2 , Y 2 }, and the incidence rates of the unobserved confounders, one may still be able to obtain useful constraints on the true treatment effect via PBA. 11 The approach is straightforward: one chooses prior probability distributions for the incidence rates of the unobserved confounder, and for each of the factors F governing the associations among the attributes, except for the treatment effect itself, F x → Y2 . One then performs iterations of drawing a set of values from these distributions. For each set, one computes the value of F x → Y2 needed to realize the observed rates In this way, one obtains a posterior distribution of possible values of the true treatment effect. The number of iterations should be sufficiently large that increasing it further does not appreciably change the shape of the posterior distribution. Here, we implement all of the above in R and perform 50 000 iterations for each scenario.

APPLICATION TO AN OBSERVATIONAL STUDY
We apply our method to a study of HD influenza vaccine effectiveness, relative to SD vaccine, against influenza and influenza-associated outcomes that was conducted within the VA patient population by Young-Xu et al. 9 In this study, a difference in the rate of hospitalization for pneumonia and influenza (P&I; ICD-9 codes 480-488) in the HD and SD arms was found in the baseline period even after matching on patient comorbidities, suggesting residual confounding by indication. This was addressed through use of PERR adjustment. We will apply the above PBA methodology to assess the robustness of this approach in this particular study.
For the baseline period, the study reported event rates (in units of events per 10 000 person-weeks) of R HD,base = 3.24 and R SD,base = 2.1. For the observation period, the rates were R HD,obs = 3.45, R SD,obs = 2.98. Thus, the RRs of HD versus SD arms in the baseline and observation periods were RR base = 1.54 and RR obs = 1.16, respectively. This made the unadjusted relative effectiveness rVE HD,obs = 1 − RR obs = − 16%, suggesting (contrary to evidence from its RCT 12 ) that the HD vaccine was less effective than SD. However, the fact that the RR in the baseline period differed significantly from 1 (RR base = 1.54) suggested confounding by indication. The PERR estimate of the RR was for an HD to SD relative effectiveness where the confidence interval was obtained using the assumption of Poisson-distributed events. We apply the causal diagram of Figure 1 to this study, with the following interpretation: starting in the baseline period, an unobserved confounder C that is present in part of the patient population-we can think of it as frailty not captured in the patients' medical records-potentially causes one or more of the following: 1. an elevated likelihood of hospitalization for pneumonia/influenza during the baseline period; 2. a modified likelihood of receiving HD rather than SD vaccine via confounding by indication: the presence of C 1 increases the likelihood that the healthcare provider (HCP) will identify the patient as being at elevated risk, and this may then affect the HCP's decision whether to prescribe SD or HD; 3. an elevated likelihood that the patient will also possess the confounder in the observation period, ie, C 2 = c 2 (a frail patient is likely to remain frail); 4. if the confounder does carry over into the observation period, then an elevated likelihood of hospitalization for pneumonia/influenza during the observation period; and 5. a reduced immune response to vaccination, resulting in a reduction in the protective effect derived from both HD and SD.
Point 5 requires further explanation. It has been shown that frailty can substantially reduce the effectiveness of influenza vaccine against influenza-associated hospitalization. 13 Furthermore, it has been shown that the relative efficacy of HD versus SD does not vary significantly between frail and nonfrail individuals. 14 We thus make the assumption that frail individuals receive a weakened vaccine effect, such that the RR due to vaccination is multiplied by a factor f > 1 for both HD and SD vaccination. It can be shown (see Appendix A) that this is equivalent to increasing F c2 → Y2 by the same factor, which, in turn, suggests that F c2 → Y2 > F c1a → Y1 . In our PBA in the following (and in some of the additional analyses in Appendix B), we use the parameterization F c → Y2 = fF c → Y1 . From Section 4, we know that, when f ≥ 1, this puts us in the regime of a pessimistic PERR, ie, one that will underestimate the relative effectiveness of HD, all other things being equal.
In the context of this study, the edge directed from Y 1 to X in Figure 1 represents the possibility that the HCP's choice of vaccine is directly influenced by whether or not a patient was hospitalized for pneumonia/influenza during the baseline period. P&I hospitalization can also influence vaccine choice by the causal path from Y 1 to C 1b to X: hospitalization may cause frailty (instead of/in addition to the other way around), and the presence of this newly acquired frailty may then influence vaccine choice. As we saw in Section 4, both the direct and indirect path can cause the PERR to be overoptimistic.
To investigate the possibility of a direct link from Y 1 to X, we conducted an interview study among nurse infection preventionists in 25 Veterans Health Administration (VHA) facilities to gain understanding of the decision-making process governing the selection of HD versus SD for a given patient within VA facilities. The study is described in Appendix C. In no case was previous hospitalization for pneumonia/influenza reported as a criterion for preferentially administering HD. On the contrary, in two (8%) of the facilities, recent hospitalization for pneumonia was reported as a possible contra-indication for HD administration. This suggests that F y1 → X ≤ 1. We conduct our probabilistic bias analyses by performing Monte Carlo simulations consisting of 50 000 realizations each (some of these are discarded due to having infeasible combinations of inputs). Each realization proceeds via a "draw and adjust" procedure, as follows.
1) A set of observed values for R( 1 |x), R( 1 |x), R( 2 |x), and R( 1 |x) is drawn from the results of the study of Young-Xu et al, assuming, as they do, that event counts are Poisson-distributed. 2) All DAG parameters apart from F c1 → X , r x , F x → Y2 and r y2 are drawn from their respective prior distributions, and any required relationships among parameters are imposed.
We choose uniform distributions for all input parameters. In Appendix B, we examine the effect of applying progressively more constraints on the input parameter ranges and relationships among them (Table B1). Corresponding posterior distributions of rVE HD,true are shown in Figure B1. As can be seen, under minimal assumptions (Analysis 1), PERR is very likely to significantly overestimate rVE HD . The distribution of rVE HD,true becomes tighter and/or moves further to the right in each successive analysis, and for Analyses 6 and 7, PERR is more likely to underestimate than overestimate the true effect size.
We now seek to identify specific constraints appropriate for this study. Table 1 gives a set of constraints, informed by published literature, on the relationship between hospitalization and frailty state transitions (see the work of Gill et al 15 ) and on hospitalization rates, including for P&I, of high-risk and low-risk individuals (see the work of Mullooly et al 16 ). In using the latter, we assume that the RR of P&I hospitalization of high-risk versus low-risk subjects can be used as a proxy for that of frail versus nonfrail subjects. The results of the PBA are given in Table 2, shown graphically in Figure 8, and compared to the PERR estimate in the study of Young-Xu et al. 9 Both distributions have a similar lower 95% CI bound, but both the median and upper bound of the PBA posterior distribution are higher than those of the PERR estimate. In Table 1, we compare the distributions using two different metrics: the common-language effect size (CLES) 17 and the two-sample Hodges-Lehmann estimator 18 for the difference between two populations. The CLES gives the probability that a sample drawn from the distribution of rVE HD,PERR is greater than a sample drawn from the distribution of rVE HD,true . In other words, it gives the probability that PERR overestimates the true effect size. The two-sample Hodges-Lehmann estimator is a nonparametric estimator of the median difference between a pair of samples drawn from the two distributions. Both comparisons suggest that the PERR estimate in the study of Young-Xu et al is more likely to have underestimated than overestimated the true relative effectiveness of HD vaccine.

DISCUSSION
In multiple studies, 1-3 PERR adjustment has been demonstrated to perform well in controlling for bias due to unmeasured confounding. The method is also appealing due to its relative simplicity, both conceptually and in implementation.   However, it has been shown that this method may in some cases increase rather than decrease bias. 5 We further explored PERR performance and used Bayesian network calculations applied to the causal diagram representation of a previously published observational study of the relative effectiveness of HD versus SD influenza vaccine in the US VA patient population. Using PBA, we showed that, in applying the PERR estimator to control for unmeasured confounding, this particular study is more likely to have underestimated rather than overestimated the true effect size. This is not to argue that the PERR should be categorically discarded in favor of PBA on Bayesian networks. Under appropriate conditions, the PERR alone will suffice. With reference to the causal diagram of Figure 1, the PERR can safely be used if all of the following apply: (i) A direct causal association between baseline period event Y 1 and treatment X can be ruled out. Such an association can masquerade as an unmeasured confounder, and naïve application of the PERR in this case acts to increase rather than decrease bias. (ii) There is no variation in the strength of the confounder effect between baseline and observation period. (iii) Individuals can neither gain nor lose the confounder during the study period. (iv) A bidirectional relationship between confounder and baseline-period event can be ruled out; that is, presence of a baseline event does not affect the probability that the confounder is present. (v) No direct causal association between baseline period event Y 1 and observation period event Y 2 exists. In practice, this is the least critical constraint, since the effect of such an association on the accuracy of the PERR estimator is modest and bounded.
If (i) or (ii) is violated, but one knows the directionality of the relationship, then, depending on the study question, the PERR may still be able to provide a useful constraint on the treatment effect. This work is subject to a number of limitations. To begin with, although the causal diagram we use is relatively generic, it by no means constitutes an adequate description of every possible type of real-world study. For example, loss of study subjects due to mortality or dropout is not taken into account. In our example of the study of Young-Xu et al, 9 no subjects were lost during the baseline period, and less than 3% were lost during the observation period (unpublished data). Uddin et al 5 did consider loss of subjects and showed that this will cause the PERR to overestimate the treatment effect. For a loss rate of less than 5%, they showed underestimation by well below 1%. However, any study in which mortality/dropout plays a more significant role will need to add this effect into the model. Also, because PERR involves repeat measurements of the same population, correlation effects may be present. In the study of Young-Xu et al, this issue, as well as the possibility of cluster effects at the VA facility level, was addressed by conducting a sensitivity analysis using generalized estimating equations. Even so, our analysis abstracts time dependence into simply a baseline and an observation period. Studies in which time evolution of subjects is more complex will need to utilize a correspondingly more complex DAG, for example, one in which baseline and observation are broken down into multiple subperiods. It should also be emphasized that, in our PBA re-analysis of the study of Young-Xu et al, the sources we use to constrain PBA inputs are derived from a targeted literature search, not an exhaustive literature review. More broadly, any analysis of this type can only be as valid as the DAG that informs it. The construction of a DAG ultimately relies on expert opinion, and no amount of expert opinion can guarantee that a causal link will not be misspecified or overlooked.
Despite the above limitations, our study offers an approach to understanding under what scenarios the PERR is likely to provide an unbiased estimate of the treatment, and a methodology to bound the bias when it is present.

DISCLOSURES
E. Thommes, R. van Aalst, J. Lee, and A. Chit are employees of Sanofi Pasteur. S. Mahmud and Y. Young-Xu have received funding from Sanofi Pasteur in the context of this and previous studies. J. Snider is an employee of and holds equity in Precision Health Economics, which has received consulting fees from Sanofi Pasteur for this and previous studies.

DATA AVAILABILITY STATEMENT
The computer code used to generate the results of this study is included in the published article's supplementary information files.

APPENDIX A
With specific reference to the HD influenza study of Young-Xu et al, 9 let us rewrite Equation (9) in terms of the absolute (rather than relative) HD and SD effect and make things explicit by redefining the domain of X as dom(X) = [HD, SD]: Now, suppose further that individuals who are frail in the observation period (C 2 = c 2 ) receive an attenuated vaccine effect, such that F HD → Y2 and F SD → Y2 are both multiplied by a factor f > 1. Thus, we have which can be rewritten to absorb f into the term involving F c → Y2 Finally, if we define r ′ 2 = r 2 F SD→Y 2 , we can write this as If, in the absence of the attenuation, F c2 → Y2 = F c1a → Y1 , then with attenuation factor f > 1,

APPENDIX B
We begin with minimal assumptions about the input parameters (Analysis 1) and then consider a series of scenarios wherein we apply progressively more restrictive constraints on the prior ranges of the input parameters. All ranges are

FIGURE B1
Results of probabilistic bias analysis to assess the true relative effectiveness of high-dose (HD) versus standard-dose (SD) influenza vaccine, rVE HD , for the scenarios in Table B1 [Colour figure can be viewed at wileyonlinelibrary.com] taken as uniform, and in all scenarios, we take r c1a ∈ [0.25,0.75], F c1a → Y1 ∈ [1,10], and F y1 → Y2 ∈ [1,10]. For each analysis, we report the distributions of rVE HD,true and compare it to rVE HD,PERR , the PERR estimate derived by Young-Xu et al. Specifically, we wish to quantify to what degree the PERR over/underestimates the true treatment effect. We do so using two different metrics: the CLES 17 ) and the two-sample Hodges-Lehmann estimator for the difference between two populations. 18 The CLES gives the probability that a sample drawn from the distribution of rVE HD,PERR is greater than a sample drawn from the distribution of rVE HD,true . In other words, it gives the probability that PERR overestimates the true effect size. The two-sample Hodges-Lehmann estimator is a nonparametric estimator of the median difference between a pair of samples drawn from the two distributions. The setups and results of the simulations are given in Table B1. Results are shown graphically in Figure B1.

APPENDIX C INTERVIEW SURVEY: HIGH-DOSE VACCINE ALLOCATION IN THE VETERANS HEALTH ADMINISTRATION BACKGROUND
The objective of this survey study was to understand the contributing factors to the variations in the HD influenza vaccine coverage among VHA facilities and to gain insight into the decision-making process for determining eligibility for and access to HD vaccine. Starting with minimal assumptions about the relationships among C

METHODS
and Y 2 , we add assumptions (described in the leftmost column) one by one and apply the corresponding constraints on the input parameters. Cells corresponding to tightened constraints are shown in light gray, and darker gray for the one instance of a constraint further tightened hospitalization has no effect on frailty determined by CPT codes (HD: 90662, SD: 90655-90659, Q2034-Q2039) from VHA EMRs. VHA medical centers were stratified into four groups according to HD vaccine proportion: 0%-4%, 5%-29%, 30%-69%, and 70%-100%. Initially, 10 facilities were randomly selected from each of the four HD vaccine proportion groups. For each site, we randomly selected and invited one nurse infection preventionist at a time to participate in a 10-15-minute telephone interview. If the interviewee was unable to provide all requested information at the time of the interview, follow-up interviews were performed with the individual, and/or additional facility personnel, such as other infection preventionists or pharmacists. If no response to the initial invitation(s) was received within a specified amount of time, another facility from within the same HD vaccine proportion group was selected at random, and the process repeated.

RR of direct
The content of the interview focused on the two most recent influenza seasons: 2015-2016 and 2016-2017. Twenty-two questions were organized into the following four areas.
HD vaccine ordered and used Nine (36%) facilities provided partial data for either the quantity of HD vaccine ordered or percentage used, or were able to provide the information for one season only. The ordered, used and percentage HD used varied greatly across the four groups and within each group (Table C2). Only 8 out of 20 facilities (40%) and 10 out of 19 (53%) reported utilizing 100% of the ordered HD in 2015-2016 and 2016-2017, respectively. Of the 17 facilities that reported ordered HD quantities for both seasons, 10 (59%) increased, 1 (6%) decreased, and 6 (35%) did not alter the amount ordered for the following 2016-2017 influenza season as compared to 2015-2016. Also, 40% (N = 10) increased the amount of HD ordered for the most recent season, 2017-2018, as compared to 2015-2016 and 2016-2017. Decision making process for HD vaccine acquisition Primary decision-makers for the quantity of seasonal HD supply reported were similar across the four groups as pharmacy and/or a multidisciplinary committee consisting of representatives from various specialties, most commonly Infection Control and Infectious Diseases (Table C3). In seven (28%) facilities, pharmacy was a sole decision-maker for determining the amount of HD to be ordered for a season.
Overall, 10 (40%) facilities reported a shortage or inability to meet the demand for HD vaccine, and, for some, this occurred despite reported unused HD vaccine at the end of the season (N = 3). In these instances, the leftover doses were attributed to strict adherence to the criteria for HD administration (N = 1), extra supply received from other facilities (N = 1), and for one, no explanation was offered. Of the 12 facilities that reported unused HD, the remaining 9 (75%) reported no shortage or unmet demand. Reasons for the surplus cited included staff unaware of HD availability (N = 3), strict adherence to the criteria for HD administration (N = 1), lack of demand (N = 1), low vaccination rates overall (N = 1), possibly poor documentation of administered HD (N = 1), and for two, no explanation was offered.
The cost of HD vaccine (N = 2), lack of conclusive evidence (N = 5), and unused doses of HD in the prior season (N = 2) were the reported barriers to increasing HD supply for the future seasons. A common theme among the surveyed infections preventions and pharmacists was the perception of a lack of clear guidance from the Centers for Disease Control and Prevention (CDC) or in the scientific literature regarding administration of HD vaccine. Criteria, policies, and practices for HD administration Twenty (80%) of surveyed facilities have a standing nursing order for HD administration and the use of such a practice varied from 100% for Group A to 67% for Group C (Table C4). Similarly, a majority or 18 (72%) have a facility-wide policy on HD encompassing: age alone (N = 13); a combination of age and presence of certain comorbidities (N = 4); or positive HIV status (N = 1). Ten (40%) facilities have department-specific policies or practices for: Home-Based Primary Care (N = 2), acute care (N = 3), and Community Living Center (CLC)-residing (N = 8), patients, as well as geriatric (N = 1), Infectious Diseases (N = 1), and HIV (N = 1) clinics. Provider discretion for HD administration on an individual basis despite the presence of facility-wide or department-specific policies was reported by 20 (80%) of facilities, and, furthermore, for 18 (72%), respondents reported that providers specifically target patients they determine to be at high-risk not otherwise covered in an existing policy.
Overall, 14 (56%) reported a facility-wide policy, 6 (24%) department-specific policies, and 4 (16%) a combination of both facility-wide and department-specific policies, while 1 (4%) reported relying solely on provider's discretion for administering HD in patients who are aged 65 years and older. All facilities in the highest HD proportion group (D) reported using a facility-wide policy on the basis of age alone and no department-specific policies; however, 60% (N = 3) of these reported that providers are able to use their discretion to target high-risk patients if they do not meet the age criterion.
Defining high-risk in the setting of HD vaccination varied from one facility to another: 22 (88%) consider advanced age (defined as 65+ (N = 20), 80+ (N = 1) or 85+ (N = 1)) to be a criterion, which aligns with the policies described above (Table C5). In 18 (72%) of these 22 facilities, meeting the age criterion alone was sufficient to be considered at high-risk, whereas for four (16%), such a consideration necessitated a combination of advanced age and the presence of certain chronic conditions such as diabetes, chronic obstructive pulmonary disease, chronic heart disease, or immunocompromised status. In addition to these, other patient criteria reported included, but were not limited to positive HIV status, CLC residency or being seen in a high-risk clinic. Contraindications for HD administration Eleven (44%) of the surveyed facilities reported contra-indications for HD administration (Table C6). Two (8%) reported that a recent pneumonia hospitalization could be a reason to give a SD instead of HD vaccine to an individual. The most commonly reported reason for delay or avoidance of HD vaccination was prior history of allergic reaction to an influenza vaccine or its components (N = 8). Presence of fever was cited as the second most common contra-indication (N = 5), followed by already immunized (N = 2), and history of Guillain-Barre syndrome (N = 2).