Skip to main content
Download PDF

Bias amplification in the g-computation algorithm for time-varying treatments: a case study of industry payments and prescription of opioid products

BMC Medical Research Methodology volume 22, Article number: 120 (2022) Cite this article

Abstract

Background

It is often challenging to determine which variables need to be included in the g-computation algorithm under the time-varying setting. Conditioning on instrumental variables (IVs) is known to introduce greater bias when there is unmeasured confounding in the point-treatment settings, and this is also true for near-IVs which are weakly associated with the outcome not through the treatment. However, it is unknown whether adjusting for (near-)IVs amplifies bias in the g-computation algorithm estimators for time-varying treatments compared to the estimators ignoring such variables. We thus aimed to compare the magnitude of bias by adjusting for (near-)IVs across their different relationships with treatments in the time-varying settings.

Methods

After showing a case study of the association between the receipt of industry payments and physicians’ opioid prescribing rate in the US, we demonstrated Monte Carlo simulation to investigate the extent to which the bias due to unmeasured confounders is amplified by adjusting for (near-)IV across several g-computation algorithms.

Results

In our simulation study, adjusting for a perfect IV of time-varying treatments in the g-computation algorithm increased bias due to unmeasured confounding, particularly when the IV had a strong relationship with the treatment. We also found the increase in bias even adjusting for near-IV when such variable had a very weak association with unmeasured confounders between the treatment and the outcome compared to its association with the time-varying treatments. Instead, this bias amplifying feature was not observed (i.e., bias due to unmeasured confounders decreased) by adjusting for near-IV when it had a stronger association with the unmeasured confounders (≥0.1 correlation coefficient in our multivariate normal setting).

Conclusion

It would be recommended to avoid adjusting for perfect IV in the g-computation algorithm to obtain a less biased estimate of the time-varying treatment effect. On the other hand, it may be recommended to include near-IV in the algorithm unless their association with unmeasured confounders is very weak. These findings would help researchers to consider the magnitude of bias when adjusting for (near-)IVs and select variables in the g-computation algorithm for the time-varying setting when they are aware of the presence of unmeasured confounding.

Peer Review reports

Background

In most epidemiologic studies, there is a dilemma between including only true confounders or all possible confounders. It is often believed that including as many covariates as possible would reduce the bias due to unmeasured confounding between the treatment and outcome. However, as suggested in prior literature, this is not always true, and conditioning on variables that are associated with the outcome only through the treatment, called instrumental variables (IVs), can increase the bias of treatment effect estimates in the point-treatment settings under the presence of unmeasured confounding [1,2,3,4,5]. This is also the case for near-IVs, variables that are weakly associated with the outcome not through the treatment, and these covariates are often called ‘bias amplifiers’ [1, 2]. For example, a previous study demonstrated that adjusting for glaucoma diagnosis (by including it in the propensity score model)—a near-IV of statin use vs. glaucoma drugs use—moved the estimated effect of statin (vs. glaucoma drug) on mortality or hip fracture risk away from the expected effect based on the results from randomized controlled trials [6]. Some previous studies have emphasized the practical challenges to determine which variables are confounders or IVs [2, 3], suggesting the need for careful consideration of including strong predictors of the treatment in the model. If the variables are weakly associated with the outcome not through the treatment (i.e., near-IVs), it might be better to present the findings from both models with and without such variables in the point-treatment settings. However, it has been still unclear whether and the extent to which this bias amplifying feature of IVs and near-IVs can be applied to the time-varying treatment settings.

G-computation, which is the computational algorithm of g-formula, is one of the methods to estimate the causal effect of time-varying treatments accounting for time-varying confounders that are affected by the treatment [7,8,9]. In practice, we need the fits of 1) the regression model of the outcome on the time-varying treatments and time-varying covariates and 2) the regression models of the time-varying covariates on previous treatments and covariates, as well as the Monte Carlo integration from these model fits in the g-computation algorithm. Note that, for time-varying treatment settings, there can be IV or near-IV for each time-point which may influence the g-computation estimator in an unpredictable way through distinct regression models in the presence of unmeasured (time-fixed or time-varying) confounders. Hence, it is crucial to quantify the influence of including IV or near-IV in the regression models on the g-computation estimates.

The goal of our study is to compare the magnitude of bias due to the adjustment of IV or near-IV across different relationships between (near-)IV and treatments in the time-varying settings. After describing a case study of industry marketing payments and prescriptions of opioid products, we conducted Monte Carlo simulation studies to investigate the extent to which the bias due to unmeasured confounders is amplified by adjusting for IV or near-IV (which could also reduce bias as a proxy of an unmeasured confounder). Our study will guide readers to consider the possible magnitude of bias due to adjustment for potential IVs, and thus help them to select covariates in their g-computation algorithms under the time-varying setting in future epidemiological research.

Methods

A case study of industry payments and prescriptions of opioid products

Opioid overdose is a major public health crisis in the United States [10]. As initial exposure to opioid prescriptions by physicians is known as a risk factor of opioid misuse and dependence [11,12,13], it is important to understand the upstream determinants of physician prescription of opioids. The financial relationship between physicians and the pharmaceutical industry has received substantial attention as it may affect physicians’ clinical practice [14, 15]. Since the launch of the Open Payments program under the Physician Payment Sunshine Act which requires the pharmaceutical industry to publicly report data on all payments and ownership interests made to licensed physicians and teaching hospitals in 2013 [16], a growing body of literature linking the financial physician-industry relationship with physicians’ prescriptions of opioids has been published [17,18,19,20,21]. All previous studies consistently showed the association between the receipt (or the number of encounters) of industry marketing payments at a single time point and an increased number of opioid prescriptions. This association was observed even among physicians who already received industry payments in the previous year [21], generating a hypothesis that physicians who received the payments at multiple time points would prescribe opioids than physicians who never received the payments. However, the evidence is lacking about the potential impact of industry marketing at multiple time points on physicians’ prescription of opioid products.

In addition, while the receipt of industry payments for non-opioids has not been adjusted for in some previous studies, it is strongly correlated with the receipt of opioid-related industry payments [21] and may also be associated with opioid prescriptions through patient-level characteristics that are not available in the Open payments data (i.e., near-IV). For example, physicians who treat patients with severe cancer may seek an educational opportunity through industry marketing of chemotherapy, and they would also have a high chance to prescribe opioids to control pain due to cancer. This possible link poses a question of whether the results change by adjusting for the receipt of industry payments for non-opioids.

Therefore, we employed a g-computation algorithm to account for the time-varying confounding due to change in opioid prescription pattern over time that was partially affected by industry marketing (i.e., physicians who received industry marketing might prescribe more opioids which would motivate the pharmaceutical industry to conduct further marketing in the following year). We then compared estimated effects between models with and without the receipt of industry payments for non-opioids, a potential (near-)IV.

Data sources and causal structure

We used data from the Centers for Medicare and Medicaid Services (CMS) Open Payments database 2015, 2016, and 2017 [22] linked with the CMS National Plan & Provider Enumeration System (NPPES) database [23], the CMS Physician Compare database [24], and the CMS Medicare Provider Utilization and Payment Data 2016, 2017, and 2018 [25]. We restricted physicians to those who treated Medicare beneficiaries and had physician-level characteristics, resulting in the final analytical sample of 250,944 physicians. Detailed information on each database can be found in previous literature [20, 21].

Throughout this paper, we let T1 and T2 denote the treatment at time 1 and time 2, respectively. We let Y denote the outcome of interest. We let X1 denote the common cause of T1, T2, and Y, let X2 denote the common cause of T2 and Y affected by T1, and let X3 denote the common cause of T1 and T2 (i.e., IV). We let Yt1,t2 denote the potential outcome if the treatment had taken values T1 = t1 and T2 = t2.

Causal diagram is shown in Fig. 1. Our treatments of interest (T1 and T2) are the receipt of general payments (all forms of non-research payment including meals, speaker compensation, honoraria, travel and lodging, consulting fees, gifts, and education materials) related to opioids in 2016 (T1) and 2017 (T2). Our outcome of interest (Y) was the opioid prescribing rate (the percentage of the total claims represented by opioid claims) in 2018. Covariates at baseline (X1) included physicians’ sex, years in practice, specialty (30 categories listed in Supplementary Table 1), attended medical school (top-20 U.S. medical schools, U.S. medical schools ranked between 21 and 50, or others based on the U.S. News & World Report research ranking), the average age of beneficiaries, the proportion of male beneficiaries, average hierarchical condition category score of beneficiaries, and opioid prescribing rate in 2016. Time-varying confounder (X2) was the opioid prescribing rate in 2017. We considered the receipt of industry payments for non-opioids in 2016 (X3) as a near-IV that was strongly associated with X1 and X2 but was weakly associated with Y only through X1 and characteristics of patients (e.g., comorbidities, socioeconomic status, etc.) who each physician treated. Because such patient-level characteristics were not available in the Open Payments data while it may influence the physicians’ receipt of industry payments in general, we evaluated whether the results changed between models with and without adjustment of the receipt of industry payments for non-opioids in 2016. More detailed information on treatment, outcome, and covariates can also be found in previous literature [20, 21].

Fig. 1
figure 1

Causal diagram assumed for example about the effect of industry marketing for opioid products on physicians’ opioid prescribing rate. Patient-level characteristics were not available in the Open Payments data

Full size image

Statistical analysis

We employed the following steps of the g-computation algorithm to estimate the mean difference in the opioid prescribing rate in 2018 according to the receipt of opioid-related industry marketing payments in 2016 and 2017. Analytical steps in this case study are shown as follows:

  1. 1)

    Fit a linear regression model to predict opioid prescribing rate in 2017 (X2) given the receipt of opioid-related payments in 2016 (T1) and baseline covariates (X1).

  2. 2)

    Fit a linear regression model to predict opioid prescribing rate in 2018 (Y) given the receipt of opioid-related payments in 2016 (T1) and 2017 (T2), baseline covariates (X1), and opioid prescribing rate in 2017 (X2).

  3. 3)

    Randomly assign the treatment status in 2016 (new T1) and 2017 (new T2) based on the proportion in the original dataset (i.e., new T1 and new T2 are marginally independent of all covariates [26]).

  4. 4)

    Predict opioid prescribing rate in 2017 (new X2) using the fitted regression model in step 1, newly assigned treatment status (new T1), and baseline covariates (X1).

  5. 5)

    Predict opioid prescribing rate in 2018 (new Y) using the fitted regression model in step 2, newly assigned treatment status (new T1 and new T2), newly predicted time-varying confounder (new X2), and baseline covariates (X1).

  6. 6)

    Estimate mean difference between distinct counterfactual marginal expectations of Y: E [Yt1, t2] – E [Y0, 0] for (t1, t2) = (1, 0), (0, 1), and (1, 1) where t1 and t2 were 1 when physicians received opioid-related payments in 2016 and 2017, respectively, and were 0 when they did not receive the payments in 2016 and 2017, respectively.

  7. 7)

    Calculate the 95% confidence interval by 200 bootstrapped samples (step 1–6).

Then, we compared the results with those when the g-computation algorithm in steps 1 and 2 additionally included the receipt of industry payments for non-opioids (X3).

Monte Carlo simulation study

Causal structure and data-generating process

We conducted two Monte Carlo simulation studies to obtain quantitative results adjusting for IV or near-IV. In these simulation studies, we aimed to evaluate the difference in biases due to unmeasured confounding between models with and without adjustment of perfect IV (scenario A) or near-IV (scenario B) of the time-varying treatments. In both scenarios, we simulated 10,000 datasets of sample size (N) = 500, 10,000, and 200,000, respectively.

In scenario A, we assumed that X3As (X3A_1, X3A_2, X3A_3, X3A_4, X3A_5, X3A_6, X3A_7, X3A_8, and X3A_9) are associated with Y only through T1 or T2, and therefore, perfect IVs (Fig. 2A). Each X3As has a different magnitude of association with T1 and T2 as shown in Table 1. We first generated X1 and X3As for subject i (=1, …, N). X1 was drawn from independent Bernoulli distributions with parameter 0.5. X3As were drawn from independent standard normal distributions; N (0, 1). Under the guidance of causal structure in Fig. 2A, we generated the status of T1, X2, T2, and Y for subject i (=1, …, N) using the following equations and values of each parameter in Table 1:

$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,t1}\right)=-10+\log (5){X}_1+{\alpha}_{3A\_1}{X}_{3A\_1}+{\alpha}_{3A\_2}{X}_{3A\_2}+{\alpha}_{3A\_3}{X}_{3A\_3}+{\alpha}_{3A\_4}{X}_{3A\_4}+{\alpha}_{3A\_5}{X}_{3A\_5}+{\alpha}_{3A\_6}{X}_{3A\_6}+{\alpha}_{3A\_7}{X}_{3A\_7}+{\alpha}_{3A\_8}{X}_{3A\_8}+{\alpha}_{3A\_9}{X}_{3A\_9}\\ {}T1i\sim \mathrm{Bernoulli}\ \left({p}_{i,t1}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,x2}\right)=-1+\log (5){T}_{1i}+\log (5){X}_1\\ {}X{2}_i\sim \mathrm{Bernoulli}\ \left({p}_{i,x2}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,t2}\right)=-10+\log (1.2){T}_1+\log (5){X}_1+\log (2){X}_{2i}+{\gamma}_{3A\_1}{X}_{3A\_1}+{\gamma}_{3A\_2}{X}_{3A\_2}+{\gamma}_{3A\_3}{X}_{3A\_3}+{\gamma}_{3A\_4}{X}_{3A\_4}+{\gamma}_{3A\_5}{X}_{3A\_5}+{\gamma}_{3A\_6}{X}_{3A\_6}+{\gamma}_{3A\_7}{X}_{3A\_7}+{\gamma}_{3A\_8}{X}_{3A\_8}+{\gamma}_{3A\_9}{X}_{3A\_9}\\ {}{T}_{2i}\sim \mathrm{Bernoulli}\ \left({p}_{i,t2}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,y}\right)=-1+\log (1.2){T}_{1i}+\log (1.5){T}_{2i}+\log (1.5){T}_{1i}{T}_{2i}+\log (5){X}_1+\log (2){X}_{2i}\ \\ {}{Y}_i\sim \mathrm{Bernoulli}\ \left({p}_{i,y}\right)\end{array}}$$
Table 1 Assigned values of parameters in the data-generating process of scenario A
Full size table
Fig. 2
figure 2

Causal diagrams for simulation studies. T1: treatment at time 1; T2: treatment at time 2; X1: Common cause of treatment at time 1, treatment at time 2, and outcome; X2: Time-varying confounder (i.e., confounder between treatment at time 2 and outcome affected by treatment at time 1); X3: Common cause of treatment at time 1 and treatment at time 2 (X3As, IV in scenario A; and X3B, near-IV in scenario B)

Full size image

In scenario B, we assumed that X3B is associated with X1 in addition to T1 and T2 (Fig. 2B). In this causal structure, as X3B is associated with Y through X1, X3B is not a perfect instrument (i.e., near-IV or proxy confounder). X1 was drawn from independent Bernoulli distributions with parameter 0.5. Under the guidance of causal structure in Fig. 2B, we generated the status of T1, X2, T2, and Y for subject (i=1,…,N) using the following equations. To detect the possible impact of adjusting for near-IV, we assumed that X3B is strongly associated with T1 and T2 (OR = 10.0) and varied its association with X1 as follows; β1 = 0.01, 0.05, 0.1, 0.2, and 0.3 and 0.3.

$${X}_{3 Bi}\sim N\ \left({\beta}_1{X}_1,1\right)$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,t1}\right)=-6+\log (2){X}_1+\log (10){X}_{3 Bi}\\ {}T{1}_i\sim \mathrm{Bernoulli}\ \left({p}_{i,t1}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,x2}\right)=-1+\log (5){T}_{1i}+\log (2){X}_1\\ {}X{2}_i\sim \mathrm{Bernoulli}\ \left({p}_{i,x2}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,t2}\right)=-7+\log (1.2){T}_{1i}+\log (5){X}_1+\log (2){X}_{2i}+\log (10){X}_{3 Bi}\\ {}{T}_{2i}\sim \mathrm{Bernoulli}\ \left({p}_{i,t2}\right)\end{array}}$$
$${\displaystyle \begin{array}{c}\mathrm{logit}\left({p}_{i,y}\right)=-0.5+\log (1.2){T}_{1i}+\log (1.5){T}_{2i}+\log (1.5){T}_{1i}{T}_{2i}+\log (2){X}_{1i}+\log (2){X}_{2i}\ \\ {}{Y}_i\sim \mathrm{Bernoulli}\ \left({p}_{i,y}\right)\end{array}}$$

Statistical analysis

Using the g-computation algorithm, we estimated mean difference between distinct counterfactual marginal expectations of Y:

$$\mathrm{E}\left[{\mathrm{Y}}^{\mathrm{t}1,\mathrm{t}2}\right]-\mathrm{E}\left[{\mathrm{Y}}^{0,0}\right]$$

for (t1, t2) = (1, 0), (0, 1), and (1, 1). We first fit two logistic regression models; 1) a model to predict X2 given T1 and an IV (one of X3As in scenario A) or a near-IV (X3B in scenario B) and 2) a model to predict Y given T1, T2, X2, and the same covariate used in the first regression model. Next, we used the regression coefficients obtained from these models to predict the values of the potential X2 and subsequently of the potential Y under a hypothetical intervention on T1 and T2.

True values of marginal expectations were approximately obtained in large (N = 10,000,000) sample generated with (Y1, 1, Y1, 0, Y0, 1, Y0, 0). We evaluated the biases of g-computation estimates across the 10,000 datasets. Then, we compared them with biases obtained in the model without adjusting for IV (in scenario A) or near-IV (in scenario B) under the presence of unmeasured confounding.

Results

A case study of industry payments and prescriptions of opioid products

Among 250,944 physicians included in this study, 10,826 (4.3%) physicians received opioid-related industry payments in 2016 only, 5773 (2.3%) physicians received opioid-related industry payments in 2017 only, and 12,558 (5.0%) physicians received opioid-related industry payments in both 2016 and 2017. Physicians’ demographic characteristics are shown in Supplementary Table 1. In Model 1 (without adjusting for the receipt of industry payments for non-opioids, X3), physicians who received industry payments for opioids in either 2016 or 2017 had a higher opioid prescribing rate in 2018 than physicians who did not receive the payments in 2016 and 2017 (physicians who received payments only in 2016, + 11.32% [95% CI, 9.87 to 12.78]; those who received payments only in 2017, + 7.43% [95% CI, 5.51 to 9.35]; and those who received payments in both 2016 and 2017, + 15.96% [95% CI, 15.24 to 16.69]; Table 2). We also found the association between the receipt of industry payments and the increased opioid prescribing rate in Model 2 (with adjusting for X3), but the estimated effect was smaller than those in Model 1 (physicians who received payments only in 2016, + 7.21% [95% CI, 3.95 to 10.47]; those who received payments only in 2017, + 3.63% [95% CI, 1.54 to 5.71]; and those who received payments in both 2016 and 2017, + 13.47% [95% CI, 12.20 to 14.73]).

Table 2 Estimated effect of industry marketing for opioid products on physicians’ opioid prescribing rate using g-computation model adjusting for time-varying confounders
Full size table

Monte Carlo simulation study

In scenario A, the approximate true mean difference between distinct counterfactual marginal expectation of Y under (t1, t2) = (1, 0), (0, 1), and (1, 1) were 8.4, 8.0, and 6.4 percentage point. Respectively. When we did not include IV in the model, the biases due to the presence of unmeasured confounding were 5.5 for E [Y1, 0] – E [Y0, 0], 4.3 for E [Y0, 1] – E [Y0, 0], and − 4.0 for E [Y1, 1] – E [Y0, 0]. The biases generally increased (i.e., away from the true value) when additionally adjusting for IV (Fig. 3, Supplementary Table 2). For example, we found the largest bias for E [Y1, 0] – E [Y0, 0] when the g-computation algorithms included IV which was strongly associated with T1. Likewise, we found the largest bias for E [Y0, 1] – E [Y0, 0] when the g-computation algorithms included IV which was strongly associated with T2. For E [Y1, 1] – E [Y0, 0], we found the largest bias when the g-computation algorithms included IV which was strongly associated with both T1 and T2. These trends were less clear when the sample size was 500 compared to when the sample size was 10,000 or 200,000. We did not find a clear trend in standard error in the models adjusting for IV in the g-computation algorithm (Supplementary Table 2).

Fig. 3
figure 3

Scenario A: Increase in biases of treatment effects in the model adjusting for IV (X3As) compared to those in the model without adjusting for IV under the presence of unmeasured confounder (X1). Bias of treatment effects was calculated by subtracting true values of marginal expectations obtained in a large (N = 10,000,000) sample from g-computation estimates across the 10,000 datasets in each situation. Each panel shows the magnitude of increase in biases when the g-computation algorithm additionally included IV (X3As) which is associated with T1 and T2 (odds ratio = 2.0, 5.0, or 10.0) under the presence of an unmeasured confounder (X1). Biases of [Y1, 1] – E [Y0, 0] (the right column) were multiplied by − 1 to provide intuitive information on the gap from the true estimates because the true estimates of [Y1, 1] – E [Y0, 0] were negative

Full size image

In scenario B, we found the bias amplification for all estimates by adjusting for near-IV (X3B) in the g-computation algorithm when X3B had a very weak association with unmeasured confounder X1 (β1 = 0.01 or 0.05) compared with its association with T1 and T2 (OR = 10) (Fig. 4, Supplementary Table 3). However, when the association between X3B and X1 was larger (β1 ≥ 0.1), we rather found the bias reduction by adjusting for X3B (as a proxy of an unmeasured confounder X1). This trend was observed in all estimates; E [Y1, 0] – E [Y0, 0], E [Y0, 1] – E [Y0, 0], and E [Y1, 1] – E [Y0, 0].

Fig. 4
figure 4

Scenario B: Comparison of bias of treatment effects between models with and without adjusting for near-IV (X3B) varying its relationship with unmeasured confounder (X1). Beta (β1 in X3Bi ~ N (β1X1, 1)) ranged from 0.01 to 0.3. Y-axis shows bias which was calculated by subtracting true values of marginal expectations obtained in a large (N = 10,000,000) sample from g-computation estimates across the 10,000 datasets in each situation. Biases of [Y1, 1] – E [Y0, 0] were multiplied by − 1 to provide intuitive information on the gap from the true estimates because the true estimates of [Y1, 1] – E [Y0, 0] were negative

Full size image

Discussion

Our simulation study showed that adjusting for IV of time-varying treatments in the g-computation algorithm increased bias due to unmeasured confounding, particularly when the IV has a strong relationship with the treatment of interest (i.e., T1 for E [Y1, 0] – E [Y0, 0], T2 for E [Y0, 1] – E [Y0, 0], and both T1 and T2 for E [Y1, 1] – E [Y0, 0]). Of note, we found the increase in bias even adjusting for near-IV when such variable had a very weak association with unmeasured confounders between the treatment and the outcome. Instead, the bias decreased after adjusting for near-IV when it had a stronger association with the unmeasured confounders (β1 ≥ 0.1, or having ≥0.1 correlation coefficient in our multivariate normal setting). Our findings were not stable with a small sample size (N = 500), indicating the importance of additional consideration for small sample bias.

Taken together, if there is a variable considered as a perfect IV for the time-varying exposure, we recommend not to include it in the model to obtain a less biased estimate. However, it is challenging to decide if the variable is a perfect IV or near-IV. Given that exposure is a collider between the candidate variable for IV and unmeasured exposure-outcome confounders, investigating the relationship between the candidate variable and the outcome conditional on the exposure would not overcome this challenging issue because a perfect IV can also be associated with the outcome through opened collider path due to conditioning on the exposure. Therefore, we need subject matter knowledge to determine whether the variable is a perfect IV or near-IV. If the candidate variable is likely to be near-IV for time-varying exposures, we would recommend including the variable in the g-computation algorithm because possible residual confounding would largely outweigh bias amplification due to including such a variable in the model, unless it has a very strong association with the exposure (e.g., OR > 10 with a unit standard-deviation increase) and a very weak association with the outcome (e.g., a difference of 5% standard-deviation between the presence and the absence of an unmeasured confounder). This net bias reduction is expected to be larger when the relationship between the near-IV and the unmeasured confounders gets stronger. Conducting analysis both with and without adjusting for the variable would also be an option to transparently convey the results [27, 28]. Because our simulation results were obtained through a specific data-generating process based on a case study, future studies are needed to validate our findings in other time-varying settings (e.g., different parameter values, different scales, continuous treatments, etc.).

In our case study, consistent with prior findings [17,18,19,20,21], we found that opioid-related industry marketing payments to physicians were associated with a higher rate of prescribing opioids in clinical practice. Moreover, using g-computation algorithms, we found that the receipt of opioid-related industry payments in 2016 was associated with opioid prescribing rate in 2018 regardless of whether they received opioid-related industry payments in 2017 or not. Although there are multiple factors that may influence opioid prescriptions such as insurance coverage, State laws, and the advent of abuse-deterrent opioid analgesics, our findings would extend the ongoing discussion about how industry marketing affects clinical practice to the time-varying setting (i.e., the possible ‘legacy effect’ of the industry payments on clinical practice) for not only opioids [17,18,19,20,21] but also other drugs such as cardiovascular drugs [29, 30] and insulin [31]. Our case study was based on data of licensed physicians in the US and did not include other professionals such as nurse practitioners who might have prescribed opioids during the study period. In addition, because we focused on the influence of industry marketing on overall opioid prescribing rates among physicians, whether the relationship varies across health practitioners and opioid types should be the subject of future research.

We also found that the estimated effects were generally larger in Model 1 (model without adjustment for the receipt of industry payments for non-opioids) than Model 2 (model with adjustment for the receipt of industry payments for non-opioids). If we assume there is no relationship between unmeasured confounders (e.g., patient-level characteristics, State laws, the development of abuse-deterrent opioid analgesics over the study period, etc.) and the receipt of industry payments for non-opioids, Model 1 would be preferred to Model 2. However, given that the association between the receipt of industry payments for non-opioids and the above-mentioned unmeasured confounders may not be very weak in reality, we assume that the estimated effects in Model 2 would be less biased estimates than those in Model 1 based on our simulation study.

Conclusions

In summary, we found that adjusting for a perfect IV may amplify the bias even under the setting of time-varying treatments. This was also the case for near-IVs only when the magnitude of its association with unmeasured confounders is much weaker than that with the time-varying treatments. These findings would help researchers to consider the magnitude of bias when adjusting for (near-)IVs and select variables in the g-computation algorithm for the time-varying setting when they are aware of the presence of unmeasured confounding.

Availability of data and materials

All data are publicly available at https://www.cms.gov.

References

  1. Pearl J. On a class of Bias-amplifying variables that endanger effect estimates. In proc 26th Conf Uncert Artif Intel (UAI 2010), P Grunwald & P Spirtes, eds Corvallis: Association for Uncetainty in Artificial Intelligence, pp 425–32. :8.

  2. Myers JA, Rassen JA, Gagne JJ, Huybrechts KF, Schneeweiss S, Rothman KJ, et al. Effects of adjusting for instrumental variables on Bias and precision of effect estimates. Am J Epidemiol. 2011;174:1213–22.

    Article  Google Scholar 

  3. Ding P, Vanderweele TJ, Robins JM. Instrumental variables as bias amplifiers with general outcome and confounding. Biometrika. 2017;104:291–302.

    Article  CAS  Google Scholar 

  4. Pearl J. Invited commentary: understanding Bias amplification. Am J Epidemiol. 2011;174:1223–7.

    Article  Google Scholar 

  5. Wooldridge JM. Should instrumental variables be used as matching variables? Res Econ. 2016;70:232–7.

    Article  Google Scholar 

  6. Patrick AR, Schneeweiss S, Brookhart MA, Glynn RJ, Rothman KJ, Avorn J, et al. The implications of propensity score variable selection strategies in pharmacoepidemiology: an empirical illustration. Pharmacoepidemiol Drug Saf. 2011;20:551–9.

    Article  Google Scholar 

  7. Robins J. A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect. Mathematical Modelling. 1986;7:1393–512.

    Article  Google Scholar 

  8. Snowden JM, Rose S, Mortimer KM. Implementation of G-computation on a simulated data set: demonstration of a causal inference technique. Am J Epidemiol. 2011;173:731–8.

    Article  Google Scholar 

  9. Daniel RM, Cousens SN, Stavola BLD, Kenward MG, Sterne J, a. C. Methods for dealing with time-dependent confounding. Stat Med. 2013;32:1584–618.

    Article  CAS  Google Scholar 

  10. Seth P, Scholl L, Rudd RA, Bacon S. Overdose deaths involving opioids, cocaine, and psychostimulants - United States, 2015-2016. MMWR Morb Mortal Wkly Rep. 2018;67:349–58.

    Article  Google Scholar 

  11. Wisniewski AM, Purdy CH, Blondell RD. The epidemiologic association between opioid prescribing, non-medical use, and emergency department visits. J Addict Dis. 2008;27:1–11.

    Article  Google Scholar 

  12. Barnett ML, Olenski AR, Jena AB. Opioid-prescribing patterns of emergency physicians and risk of long-term use. N Engl J Med. 2017;376:663–73.

    Article  Google Scholar 

  13. Brat GA, Agniel D, Beam A, Yorkgitis B, Bicket M, Homer M, et al. Postsurgical prescriptions for opioid naive patients and association with overdose and misuse: retrospective cohort study. BMJ. 2018;360:j5790.

    Article  Google Scholar 

  14. Spurling GK, Mansfield PR, Montgomery BD, Lexchin J, Doust J, Othman N, et al. Information from pharmaceutical companies and the quality, quantity, and cost of physicians’ prescribing: a systematic review. PLOS medicine. Public library of Science. 2010;7:e1000352.

  15. Mitchell AP, Trivedi NU, Gennarelli RL, Chimonas S, Tabatabai SM, Goldberg J, et al. Are financial payments from the pharmaceutical industry associated with physician prescribing? Ann intern med. Am Coll Physicians. 2020;174:353–61.

    Google Scholar 

  16. Kirschner NM, Sulmasy LS, Kesselheim AS. Health policy basics: the physician payment sunshine act and the open payments program. Ann Intern Med. 2014;161:519–21.

    Article  Google Scholar 

  17. Hadland SE, Cerdá M, Li Y, Krieger MS, Marshall BDL. Association of Pharmaceutical Industry Marketing of opioid products to physicians with subsequent opioid prescribing. JAMA Intern Med. 2018;178:861–3.

    Article  Google Scholar 

  18. Zezza MA, Bachhuber MA. Payments from drug companies to physicians are associated with higher volume and more expensive opioid analgesic prescribing. PLoS One. 2018;13:e0209383.

    Article  Google Scholar 

  19. Nguyen TD, Bradford WD, Simon KI. Pharmaceutical payments to physicians may increase prescribing for opioids. Addiction. 2019;114:1051–9.

    Article  Google Scholar 

  20. Hadland SE, Rivera-Aguirre A, Marshall BDL, Cerdá M. Association of Pharmaceutical Industry Marketing of opioid products with mortality from opioid-related overdoses. JAMA Netw Open. 2019;2:e186007.

    Article  Google Scholar 

  21. Inoue K, Figueroa JF, Orav EJ, Tsugawa Y. Association between industry payments for opioid products and physicians’ prescription of opioids: observational study with propensity-score matching. J Epidemiol Community Health. 2020;74:647–54.

    PubMed  Google Scholar 

  22. Centers for Medicare & Medicaid Services. Open Payments [Internet]. [cited 2019 Dec 17]. Available from: https://www.cms.gov/OpenPayments.

  23. Centers for Medicare & Medicaid Services. NPPES NPI Registry [Internet]. [cited 2019 Dec 17]. Available from: https://npiregistry.cms.hhs.gov/.

  24. Centers for Medicare & Medicaid Services. Physician Compare datasets [Internet]. Data.Medicare.Gov. [cited 2019 Dec 17]. Available from: https://data.medicare.gov/data/physician-compare.

  25. Centers for Medicare & Medicaid Services. Medicare Provider Utilization and Payment Data: Part D Prescriber Data [Internet]. [cited 2019 Dec 17]. Available from: https://www.cms.gov/Research-Statistics-Data-and-Systems/Statistics-Trends-and-Reports/Medicare-Provider-Charge-Data/Part-D-Prescriber.

  26. Wang A, Arah OA. G-computation demonstration in causal mediation analysis. Eur J Epidemiol. 2015;30:1119–27.

    Article  Google Scholar 

  27. Pimentel SD, Small DS, Rosenbaum PR. Constructed second control groups and attenuation of unmeasured biases. J Am Stat Assoc Taylor & Francis. 2016;111:1157–67.

    Article  CAS  Google Scholar 

  28. Brookhart MA, Stürmer T, Glynn RJ, Rassen J, Schneeweiss S. Confounding control in healthcare database research: challenges and potential approaches. Med Care. 2010;48:S114–20.

    Article  Google Scholar 

  29. Inoue K, Figueroa JF, DeJong C, Tsugawa Y, Orav EJ, Shen C, et al. Association between industry marketing payments and prescriptions for PCSK9 (Proprotein convertase Subtilisin/Kexin type 9) inhibitors in the United States. Circ: Cardiovasc Qual Outcomes Am Heart Assoc. 2021;14:e007521.

  30. DeJong C, Aguilar T, Tseng C-W, Lin GA, Boscardin WJ, Dudley RA. Pharmaceutical industry-sponsored meals and physician prescribing patterns for Medicare beneficiaries. JAMA Intern Med. 2016;176:1114–22.

    Article  Google Scholar 

  31. Inoue K, Tsugawa Y, Mangione CM, Duru OK. Association between industry payments and prescriptions of long-acting insulin: an observational study with propensity score matching. PLOS medicine. Public library of. Science. 2021:18:e1003645.

Download references

Acknowledgements

Not applicable.

Funding

KI was supported by National Institutes of Health (NIH)/NIDDK grant F99DK126119, Honjo International Foundation Scholarship, the Japan Society for the Promotion of Science (JSPS; 21 K20900), and the Program for the Development of Next-generation Leading Scientists with Global Insight (L-INSIGHT), sponsored by the Ministry of Education, Culture, Sports, Science and Technology, Japan. This article does not necessarily represent the views and policies of the NIH. TS was supported by JSPS (20 K11716). The funders had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Author information

Authors and Affiliations

  1. Department of Social Epidemiology, Graduate School of Medicine, Kyoto University, Floor 2, Science Frontier Laboratory, Yoshida-konoe-cho, Sakyo-ku, Kyoto, 604-8146, Japan

    Kosuke Inoue & Naoki Kondo

  2. Department of Epidemiology, UCLA Fielding School of Public Health, Los Angeles, CA, USA

    Kosuke Inoue

  3. Department of Health Data Science, Graduate School of Data Science, Yokohama City University, Yokohama, Kanagawa, Japan

    Atsushi Goto

  4. Institute for Future Initiatives, The University of Tokyo, Tokyo, Japan

    Naoki Kondo

  5. Department of Information and Computer Technology, Faculty of Engineering, Tokyo University of Science, Tokyo, Japan

    Tomohiro Shinozaki

Authors
  1. Kosuke Inoue
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Atsushi Goto
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Naoki Kondo
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Tomohiro Shinozaki
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

KI developed the study concept, analyzed and interpreted the data, supervised, and wrote the manuscript. AG, and TS developed the study concept, interpreted the data, and reviewed/edited the manuscript. NK contributed to the discussion and reviewed/edited the manuscript. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Kosuke Inoue.

Ethics declarations

Ethics approval and consent to participate

The UCLA IRB review was not required for the present study (IRB#18–001960) because we only used publicly available data without identifiable information. All methods were performed in accordance with the Declaration of Helsinki.

Consent for publication

Not applicable as this is a secondary data analysis project.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Inoue, K., Goto, A., Kondo, N. et al. Bias amplification in the g-computation algorithm for time-varying treatments: a case study of industry payments and prescription of opioid products. BMC Med Res Methodol 22, 120 (2022). https://doi.org/10.1186/s12874-022-01563-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12874-022-01563-3

Keywords

BMC Medical Research Methodology

ISSN: 1471-2288

Contact us