A novel weighting method to remove bias from within-subject exposure dependency in case-crossover studies

Abstract

Background

Case-crossover studies have been widely used in various fields including pharmacoepidemiology. Vines and Farrington indicated in 2001 that when within-subject exposure dependency exists, conditional logistic regression can be biased. However, this bias has not been well studied.

Methods

We have extended findings by Vines and Farrington to develop a weighting method for the case-crossover study which removes bias from within-subject exposure dependency. Our method calculates the exposure probability at the case period in the case-crossover study which is used to weight the likelihood formulae presented by Greenland in 1999. We simulated data for the population with a disease where most patients receive a cyclic treatment pattern with within-subject exposure dependency but no time trends while some patients stop and start treatment. Finally, the method was applied to real-world data from Japan to study the association between celecoxib and peripheral edema and to study the association between selective serotonin reuptake inhibitor (SSRI) and hip fracture in Australia.

Results

When the simulated rate ratio of the outcome was 4.0 in a case-crossover study with no time-varying confounder, the proposed weighting method and the Mantel-Haenszel odds ratio reproduced the true rate ratio. When a time-varying confounder existed, the Mantel-Haenszel method was biased but the weighting method was not. When more than one control period was used, standard conditional logistic regression was biased either with or without time-varying confounding and the bias increased (up to 8.7) when the study period was extended. In real-world analysis with a binary exposure variable in Japan and Australia, the point estimate of the odds ratio (around 2.5 for the association between celecoxib and peripheral edema and around 1.6 between SSRI and hip fracture) by our weighting method was equal to the Mantel-Haenszel odds ratio and stable compared with standard conditional logistic regression.

Conclusion

Case-crossover studies may be biased from within-subject exposure dependency, even without exposure time trends. This bias can be identified by comparing the odds ratio by the Mantel-Haenszel method and that by standard conditional logistic regression. We recommend using our proposed method which removes bias from within-subject exposure dependency and can account for time-varying confounders.

Background

The case-crossover design has been widely used since it was proposed in 1991 [1]. The design has been used in various fields including pharmacoepidemiology [2], occupational epidemiology [3], studies on traffic safety [4] and air pollution health effects [5, 6]. In case-crossover studies, individuals who have experienced the outcome (cases) act as their own controls by including one or more periods before the onset of the outcome. The period including the outcome is the case period, while period(s) prior to the case period act as the controls. The number of control periods can be large: for example, in the original article [1], one analysis involved 8766 control periods. The effect of the exposure should be brief; exposure in any period should affect the outcome in that period only, without any ‘carryover effect’ [1, 6, 7]. In addition, in case-crossover studies, the rate of outcome occurrence is usually assumed to be unchanged during exposed or unexposed periods, respectively. However, like other case-only studies, the case-crossover study has an advantage that the effect of time-invariant confounders is automatically controlled because the case period is compared with control period(s) of the same individual [7,8,9]. The case-crossover study also has unique characteristics. For example, the original unidirectional case-crossover study does not include periods after the outcome occurs. Thus, there is no bias due to the outcome influencing future exposures or future observation periods, unlike other case-only studies such as self-controlled case series [10,11,12].

However, case-crossover studies are susceptible to at least two types of major biases. The first is bias due to time trends in the exposure which can be removed using a variant of case-crossover studies, the case-time-control design [9, 13, 14]. The second is bias due to within-subject exposure dependency or autocorrelation in an individual’s exposure history [6, 15, 16]. This bias is particularly important in pharmacoepidemiology because drug use on 1 day is rarely independent from use in the preceding days. However, the potential for this bias has had little attention in the pharmacoepidemiology literature, and standard conditional logistic regression has been used without assessing whether bias due to within-subject exposure dependency exists in many case-crossover studies [17,18,19,20,21]. This may be in part because this bias has been considered rather minor when exposure is stationary [15].

In this paper, we show that within-subject exposure dependency may produce unignorably large bias even if exposure is stationary. We also describe how to remove bias from within-subject exposure dependency, when time trends in the exposure do not exist.

Motivating examples

Our motivation for investigating bias in case-crossover studies due to within-subject exposure dependency was a hypothetical drug treatment with a specific cyclic exposure pattern, where standard conditional logistic regression may produce bias. We also outline real-world pharmacoepidemiological studies in Japan and Australia where our proposed method is used.

Drug treatment with specific exposure pattern

We illustrate our motivating example using a drug treatment with a cyclic pattern consisting of two periods of treatment followed by one period with no treatment. The exposure pattern clearly shows autocorrelation, with the probability of exposure in one treatment period dependent on previous periods. We show that estimates of the odds ratio for the exposure-outcome association may be biased, depending on the number of control periods chosen.

In a hypothetical unidirectional case-crossover study with 4 periods (one case period and 3 control periods), only 3 exposure patterns are possible. If exposure and non-exposure are indicated by 1 and 0, respectively, the 3 exposure patterns are (1101), (1011), and (0110). Therefore, subjects with an exposed case period can have only 1 unexposed and 3 exposed periods, while those with unexposed case period can have only 2 unexposed and 2 exposed periods. In such a case-crossover study, when the true odds ratio (OR) is 4.0, the OR is underestimated as 3.2 using standard conditional logistic regression for matched case-control studies.

However, if there are two control periods, the possible exposure patterns are (110), (101), and (011). In this case, irrespective of whether the case period is exposed or unexposed, all subjects have 1 unexposed and 2 exposed periods and there is no bias in the estimate of the OR. Similarly, there is no bias when the total number of periods (including the case period) is an integral number of 3 (the cyclic pattern).

Overview of real-world examples

We describe two real world pharmacoepidemiological studies which will be analyzed by our proposed method. The first is an investigation into a previously reported association between celecoxib and peripheral edema [22] using Japanese data from 25 corporate-type health insurance plans [23]. From 99,821 new users of celecoxib aged between 20 and 74 years from May 2013 to April 2018, who used celecoxib after at least 180 days of non-use, we selected 311 cases who experienced peripheral edema and had both exposed and unexposed days during an 84-day study period.

The second example is an Australian study of the association between hip fracture and psychoactive medicines [24, 25]. The data were obtained from the Australian Department of Veterans’ Affairs administrative claims database. Psychoactive medicines included benzodiazepines, selective serotonin re-uptake inhibitors (SSRIs), opioids, antipsychotics and tricyclic antidepressants. The hip fracture cases were 8828 patients aged over 65 years who were hospitalized between 2009 and 2012. A previous case-crossover study of this cohort found an increased risk of hip fracture for opioids, SSRIs and antipsychotics [24]. A related case-control study found an association between hip fracture and SSRIs when used concurrently with other psychoactive medicines [25].

Methods

To avoid bias due to within-subject exposure dependency, the likelihood for estimating the odds ratio for the exposure-outcome association may be modified. At least two approaches are effective to modify the likelihood and both involve two-step weighting procedures. In the first approach, the probability of each exposure permutation is estimated in Step 1 and used as weights in the likelihood in Step 2. In the special case where this probability is the same for every permutation (“global exchangeability”), the modified likelihood is equivalent to that for the standard conditional logistic regression.

In this paper, we propose another approach that can be used if “pairwise exchangeability” is assumed. Pairwise exchangeability holds where there are no time trends in the exposure i.e., the proportion exposed in the population is stationary over the study period [15]. In this scenario, the probabilities that the case period is exposed and unexposed are estimated in Step 1 and used as weights in Step 2. When there is pairwise exchangeability, an unbiased OR can be obtained by the Mantel-Haenszel method as well.

BIAS due to conditional logistic regression for case-crossover studies

In 2001, Vines and Farrington proposed the likelihood for case-crossover studies as [15]:

$$L=\prod_{i=1}^N\frac{\exp \left(\beta {x}_{i0}\right){\sum}_{\kappa \ast }P\left\{{X}_{i0}={x}_{i\kappa \ast (0)},---,{X}_{i M}={x}_{i\kappa \ast (M)}\right\}}{\sum_{\kappa}\exp \left(\beta {x}_{i\kappa (0)}\right)P\left\{{X}_{i0}={x}_{i\kappa (0)},---,{X}_{i M}={x}_{i\kappa (M)}\right\}}$$

(1)

where X_i0 is the exposure level at the case period (m = 0) and X_im is the exposure level at the m-th control period at t (t =  − m: m = 1, 2, −-, M), x_i0 denotes the observed exposure level at the case period and x_im denotes the observed exposure level at the m-th control period, the sum in the denominator ranges over all permutations of κ of the integers {0, 1, −--, M} and the sum in the numerator ranges over the subset of the permutations for which x_iκ(0) = x_i0 of the individual i. In Eq. (1), X_i0 may be a binary exposure variable but can be a multi-level exposure variable or a vector of an exposure and time-varying confounders. When X_im denotes a binary exposure, the denominator in Eq. (1) becomes $\exp \left(\beta \right){\sum}_{\kappa_1}P\left({X}_{i0}=1,{N}_{exposed}=k\right)+{\sum}_{\kappa_0}P\left({X}_{i0}=0,\kern0.5em {N}_{exposed}=k\right)$ where N_exposed is the total number of exposed (case and control) periods, given by ${N}_{exposed}={\sum}_{m=0}^M{X}_{im}$ and ${\sum}_{\kappa_1}P\left({X}_{i0}=1,{N}_{exposed}=k\right)$ is the sum of probabilities for all the permutations of k exposed and (M + 1-k) unexposed periods where X_i0 = 1, and ${\sum}_{\kappa_0}P\left({X}_{i0}=0,\kern0.5em {N}_{exposed}=k\right)$ is that where X_i0 = 0. Data on N_exposed = k is informative only when the positivity (non-zero probability) condition is satisfied or ∑P(X_i0 = 1, N_exposed = k) > 0 and ∑P(X_i0 = 0 N_exposed = k) > 0. Otherwise, they do not contribute to the estimation of exp(β).

Let OR_VF be the estimate of exp(β) obtained from Eq. (1). The likelihood in Eq. (1) and OR_VF are in general different from the following likelihood for the standard conditional logistic (SCL) regression for individually matched case-control studies and its estimate, OR_SCL.

$$L=\prod_{i=1}^N\frac{\exp \left(\beta {x}_{i0}\right)}{\sum_{j=0}^M\exp \left(\beta {x}_{ij}\right)}$$

(2)

Vines and Farrington showed that the likelihoods in Eqs. (1) and (2) are equivalent if P{X_i0 = x_i0, − − −, X_iM = x_iM} = P{X_i0 = x_iκ(0), − − −, X_iM = x_iK(M)} for all permutations κ of {0, 1, −---, M}, that is, global exchangeability holds. For example, Eqs. (1) and (2) are equivalent when the exposure status in one period is independent from the status in any other periods and the exposure probability is the same in all of case and control periods (Appendix 1, Additional File 1). If global exchangeability does not hold, OR_SCL can be biased.

Vines and Farrington did not show explicitly how to estimate P{X_i0 = x_iκ(0), − − −, X_iM = x_iκ(M)} in Eq. (1). These probabilities may be estimated by the proportion of each permutation in the population which contain cases, or in samples representing the population such as time-controls in the case-time-control design proposed by Suissa [14]. In the next section, however, we introduce a different approach to remove the bias from within-subject exposure dependency by assuming that pairwise exchangeability is satisfied but global exchangeability may not necessarily hold.

Weighting method to remove BIAS due to within-exposure dependency

Case-crossover studies with a binary exposure

When pairwise exchangeability, P{X_i0 = 1, X_im = 0} = P(X_i0 = 0, X_im = 1} is satisfied for a binary exposure in all control periods m (m = 1, 2, −-, M), the estimate of exp(β) using the Mantel-Haenszel method (OR_MH) is unbiased whether within-subject exposure dependency exists or not [15]. In line with this finding, we will show that when pairwise exchangeability holds, the probabilities that the individual is unexposed (π₀) and exposed (π₁) at the case period, can be estimated from the cases in a case-crossover study (without requiring data from the population or time-controls). Once π₀, π₁, and the relative exposure probability π₁₀ = π₁ / π₀ are estimated, the following likelihood for case-crossover studies proposed by Greenland [26] can be used to obtain an unbiased estimate of exp(β), defined as OR_G:

$$\mathrm{L}=\prod_{i=1}^N\frac{\exp \left({\beta x}_{ic}\right){\pi}_{ic}}{\sum_k\exp \left(\beta {x}_k\right){\pi}_{ik}}=\prod_{i=1}^N\frac{\exp \left({\beta x}_{ic}\right){\pi}_c}{\sum_k\exp \left(\beta {x}_k\right){\pi}_k}=\prod_{i=1}^N\frac{\exp \left({\beta x}_{ic}\right){\pi}_{c0}}{\sum_k\exp \left(\beta {x}_k\right){\pi}_{k0}}$$

(3)

In the left-hand side of Eq. (3), π_ik is the probability that individual i has the k-th exposure level at the case period (k = 0, 1 for binary exposure), π_ic is π_ik observed, and x_ic is the exposure level observed when individual i has the outcome. In the middle of Eq. (3), subscript i is omitted in the exposure probabilities as π_ik is replaced by the expected value in the population in the current study. In the right-hand side of Eq. (3) π_k0 = π_k/π₀. Eq. (3) stands for the model where x_ik is a binary variable, multi-level exposure variable, or a vector of the observed exposure and time-varying confounders.

For a binary exposure, the right-hand side of Eq. (3) can be rewritten as

$$\mathrm{L}=\prod_{i=1}^N\frac{\exp \left({\beta x}_{ic}\right){\pi}_{c0}}{1+\exp \left(\beta \right){\pi}_{10}}$$

(4)

where x_ic = 1 and π_c0 = π₁₀ = π₁/π₀ when the individual i was exposed at the case period and x_ic = 0 and π_c0 = π₀₀ = π₀/π₀ = 1 when unexposed. As Vines and Farrington showed [15], Greenland did not estimate π_k for case-crossover studies in Eq. (3) when within-subject exposure dependency exists.

We outline a novel weighting method using a modified version of Greenland’s likelihood. We propose that π_k can be estimated, with or without within-subject exposure dependency, by assuming pairwise exchangeability. Let P_kl[m] denote the joint probability that the subject has the k-th exposure level at the case period and has the l-th exposure level at the m-th control period:

$${\mathrm{P}}_{kl\left[m\right]}=\mathrm{P}\left({X}_0={x}_k,\kern0.5em {X}_m={x}_l\right)\kern0.50em \mathrm{m}=1,2,---,\mathrm{M}$$

(5)

In Eq. (5), X₀ is the exposure status at the case period and X_m is the exposure status at the m-th control period. When the exposure variable is binary, both X₀ and X_m have two levels (k = 0, 1 and l = 0,1) and pairwise exchangeability is equivalent to stationary exposure (no time trends): when the exposure process is stationary,  P_10[m] + P_11[m] = P_01[m] + P_11[m] and this relationship leads to the pairwise exchangeability condition P_10[m] = P_01[m] (m = 1, 2, −-, M).

Using conditional probabilities, pairwise exchangeability, P{X_i0 = 1, X_im = 0} = P(X_i0 = 0, X_im = 1} can be rewritten as:

$${\uppi}_1\mathrm{P}\left(\ {X}_m=0|{X}_0=1\right)={\uppi}_0\mathrm{P}\left(\ {X}_m=1|{X}_0=0\right)$$

(6)

where π₀ = P(X₀ = 0) and π₁ = P(X₀ = 1). When both sides of Eq. (6) are summed up over M control periods (m = 1, 2, −-, M), we obtain:

$${\uppi}_1{\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=0|{X}_0=1\right)={\uppi}_0{\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=1|{X}_0=0\right)$$

(7)

The quantity ${\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=0|{X}_0=1\right)$ can be estimated by the average number of unexposed control periods (where X_m = 0) in those exposed at the case period (X₀ = 1). Similarly, the quantity ${\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=1|{X}_0=0\right)$ can be estimated as the average number of exposed control periods (X_m = 1) in those unexposed at the case period (X₀ = 0). This average, defined as $\overline{PT_{10}}$ and $\overline{PT_{01}}$, respectively, can be written as:

$$\overline{PT_{10}}={\sum}_i{PT}_{10i}/{a}_1\kern0.5em and\kern0.5em \overline{PT_{01}}\ {\sum}_i{PT}_{01i}/{a}_0$$

(8)

where PT_10i is the number of unexposed control periods (person-time) of case i who is exposed at the case period, PT_01i is the number of exposed control periods (person-time) of case i who is unexposed at the case period, and a₁ is the number of discordant exposed cases with at least one unexposed control period and a₀ is the number of discordant unexposed cases with at least one exposed control period. When ${\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=0|{X}_0=1\right)$ and ${\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=1|{X}_0=0\right)$ in Eq. (7) are substituted by $\overline{PT_{10}}$ and $\overline{PT_{01}}$, respectively, we obtain:

$${\pi}_{10}={\pi}_1/{\pi}_0={\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=1|{X}_0=0\right)/{\sum}_{m=1}^M\mathrm{P}\left(\ {X}_m=0|{X}_0=1\right)=\overline{PT_{01}}/\overline{PT_{10}}$$

(9)

When standard statistical software is used for conditional logistic regression analysis of case-crossover studies, we introduce weighting to ensure the denominator in Eq. (2) equals that in Eq. (4). Every exposed and unexposed period in case i should be weighted by w_i1 and by w_i0, respectively, defined as follows:

$${\mathrm{w}}_{i1}={\pi}_{10}/{m}_i^1\ and\ {\mathrm{w}}_{i0}=1/{m}_i^0\kern0.75em$$

(10)

where ${m}_i^1$ is the number of exposed periods and ${m}_i^0$ is the number of the unexposed periods (including both case and control periods) in case i and π₁₀ is estimated from Eq. (9). In most statistical software, w_ik (k = 0,1) may be specified as an offset variable in conditional logistic regression which uses the following likelihood:

$$L=\prod_{i=1}^N\frac{\exp \left(\beta {x}_{i0}\right)}{\sum_{j=0}^M{\mathrm{w}}_{ij}\exp \left(\beta {x}_{ij}\right)}$$

(11)

where w_ij = w_i1 and w_ij = w_i0 when j-th period is exposed and unexposed, respectively, in case i.

Using a₁ and a₀, Eq. (4) can be rewritten as

$$\mathrm{L}={\left(\frac{\exp \left(\beta \right){\pi}_{10}}{1+{\exp}\left(\beta \right){\pi}_{10}}\right)}^{a_1}{\left(\frac{1}{1+{\exp}\left(\beta \right){\pi}_{10}}\right)}^{a_0}$$

(12)

Equation (12) gives (see Additional File 1, Appendix 2) the following maximum likelihood estimate for OR_G:

$${OR}_G=\exp \left(\beta \right)=\frac{a_1}{a_0}\frac{1}{\pi_{10}}$$

(13)

When π₁₀ estimated by Eq. (9) is considered as a constant, the variance is given by:

$$\mathrm{v}\left(\upbeta \right)=\frac{1}{a_0}+\frac{1}{a_1}$$

(14)

In order to allow for the variance of the weights estimated by Eq. (9), we recommend the use of bootstrapping to estimate the 95% confidence interval using the 2.5 to 97.5 percentiles of OR_G.

From Eqs. (8), (9) and (13) we obtain:

$${OR}_G=\exp \left(\upbeta \right)=\frac{a_1}{a_0}\frac{1}{\pi_{10}}=\frac{a_1\ \overline{PT_{10}}}{a_0\ \overline{PT_{01}}}=\frac{\sum_i{PT}_{10i}}{\sum_i{PT}_{01i}}={OR}_{MH}$$

(15)

Equation (15) shows that when the model involves only one binary exposure variable, the point estimate of OR_G is the same as OR_MH.

Case-crossover studies with a binary exposure and a binary time-varying confounder

The Mantel-Haenszel estimator is unbiased for binary exposures when there is pairwise exchangeability. When there is a binary time varying confounder (z), Maclure in his original proposal of the case-crossover design recommended further stratification when using the Mantel-Haenszel method [9]. In a simulation, we followed this recommendation and stratified each subject by z. As a result, control periods where z = 0 were excluded from subjects with z = 1 at the case period. Similarly, control periods where z = 1 were excluded from those with z = 0 at the case period. We have shown that excluding these control periods can lead to bias (see Tables 3 and 4 in Results section).

On the other hand, our method can be extended to studies with time-varying confounder to analyze data of all case and control periods and all periods can be included in the analysis. For example, when there is a binary exposure (x) and a binary time-varying confounder variable (z), βx_ij in the likelihood in Eq. (11) is specified as (β, γ)(x_ij, z_ij)^T, which is equal to 0 when (x_ij, z_ij) =(0,0), β when (x_ij, z_ij) =(1,0), γ when (x_ij, z_ij) =(0,1), and β + γ when (x_ij, z_ij) = (1,1) where exp(γ) is an estimate of the odds ratio of z. Similarly to the finding that OR_SCL in Eq. (2) is unbiased when within-subject exposure dependency does not exist (Appendix 1, Additional File 1), we may estimate an unbiased OR_G in the model involving the exposure and time-varying confounder by calculating the weight from the exposure variable x (Eq. 10), if within-subject dependency does not exist for z during exposed periods (where the probability that the confounder is positive is f₁) and during unexposed periods (f₀) and the confounder is adjusted for as in the standard conditional logistic regression. Similarly as for a case-crossover study with a binary exposure, we recommend bootstrapped 2.5 to 97.5 percentiles of OR_G to estimating the 95% CI to allow for the variance of the weights (see Appendix 3, Additional File 1 for the detail; as to the relevant SAS codes, see 4-1 h and 4-2f in Appendix 4, Additional File 1).

Simulation studies

We simulated data relevant to drug therapy with no time trends and with autocorrelated exposure patterns (within-subject dependency) (Appendix 4, Additional File 1). The simulated data is created based on the following observations (i) drug treatment sometimes has a cyclic pattern which often produces within-subject exposure dependency, (ii) some outcomes (e.g., acute adverse events) tend to occur soon after the treatment is initiated but they may also occur later during the drug therapy, and (iii) some patients stop treatment for various reasons while some patients start treatment. When the rate of stopping treatment is the same as that of starting treatment, the stationarity of the exposure may be maintained in the population as follows. We simulated scenarios with and without time varying confounding.

Assume that drug treatment involves a cyclic pattern where one cycle consists of 7 days and a patient has a drug on days 1 and 4 but no drug on days 2, 3, 5, 6, and 7. Figure 1 depicts three cycles of drug treatment with 8 subgroups consisting of 1 case period and 21 control periods, where 1 period is 1 day. Subgroup A represents stoppers who stop treatment at the case period while Subgroup H represents new users who start treatment at the case period. Subgroups B to G represent patients being treated with a different timing relative to the start of the treatment cycle. In Fig. 1, the proportion of those exposed to a drug in the population is always ¼, indicating stationarity (no time trends).

Fig. 1

Exposure pattern in a hypothetical dynamic population of 80,000 patients. Patients receiving drug treatment with a cycle of 7 days are divided into 8 Subgroups A to H where 1 period is defined as 1 day. Subgroup A represents stoppers that stop treatment, Subgroup H represents starters that start treatment, and Subgroups B to G represent different exposure patterns of those being treated at the case period. The bolded outline indicates that 21 control periods can be divided into 3 cycles with the same exposure pattern. N: size of each subgroup; c0: exposure status at the case period; cm (m = 1, 2, −--, 21): exposure status at the m-th control period; Tx: treatment

Figure 2a shows 140 cases who had the outcome at the case period when the size of each of Subgroups A to H, N = 10,000, the event rate in an unexposed period (r₀) is 0.001 per period and the rate ratio is 4. We assume that the expected number of cases is determined by exposure at the case period only, or Nr₀ when unexposed and N RR r₀ when exposed at the case period.

Fig. 2

Expected frequency of cases by subgroup in Fig. 1. (a) 140 cases who had the outcome at the case period involving one binary exposure variable only, where the incidence rate per period at unexposed period (r0) is assumed to be 0.001 and the rate ratio of the exposure (RR) is assumed to be 4. (b) 184 cases who had the outcome at the case period involving one binary exposure variable and one binary time-varying confounder where the incidence rate per period at unexposed period without confounder (r0) is assumed to be 0.001, the rate ratio of the exposure (RR) is assumed to be 4, and that of the time-varying confounder (RRz) is assumed to be 2. The proportion of the time-varying confounder at unexposed periods (f0) is 0.2 and that at exposed periods (f1) in the population 0.4. The status of the confounder at control periods is randomly assigned . n: frequency of cases belonging to each ID_Subgroup in Fig. 1; c0: exposure status at the case period; z0: status of time-varying confounder at the case period

Figure 2b shows 184 cases who had the outcome at the case period when N = 10,000, the event rate in the unexposed period without time-varying confounding is 0.001 per period and the rate ratio is 4 and 2 for the exposure and time-varying confounder, respectively. The status of the time-varying confounder in the unexposed or exposed period is related to the exposure status of that period only, and f₀ = 0.2 and f₁ = 0.4 where f₀ and f₁ are the expected values of the proportion of exposed periods and unexposed periods in the population, respectively, when the confounder is positive (See Appendix 3, Additional File 1). We assume that the number of cases is as expected, or Nr₀(1 − f₀), N RR_z r₀f₀, N RR r₀(1 − f₁), and N RR RR_z r₀f₁, when the combination of the exposure and time-varying confounder variables (x, z) at the case period is (0, 0), (0, 1), (1, 0), and (1, 1), respectively. The status of the time-varying confounder in the control periods was randomly generated and is therefore independent of the exposure or time-varying confounder at different periods.

To determine the effect of the length of each time period on the estimated odds ratio, we also analyzed the data by dividing 22 days in Fig. 1 into 11, 7, 5, 3, and 2 periods (M = 10, 6, 4, 2, and 1) where 1 period included 2, 3, 4, 7 and 11 days, respectively. The status of the exposure and time varying confounder was defined by the last day of each period (Definition I in Table 1). We analyzed the data by standard conditional regression, the Vines and Farrington method, Mantel-Haenszel methods, and our modified Greenland’s method for the fixed study period of 22 days. With time-varying confounder (z), case and control periods of each subject were further stratified by z in the Mantel-Haenszel method. We used bootstrapping to estimate 2.5 to 97.5 percentiles for OR_G. Data simulation and analyses were performed using SAS 9.4 and SAS codes including those for bootstrap method are shown in Appendix 4 in Additional File 1.

Table 1 Exposure Definitions

Case-crossover studies of real-world data

The method was also applied to data from Japanese and Australian databases. The Japanese study on the association between celecoxib and peripheral edema from a previous study [22] was approved by the ethics committee of Tokyo University of Science (approval number: 18023) where obtaining the informed consent from study subjects was waived for the current study. The Japanese data came from 25 corporate-type health insurance plans extracted from Cross-Fact database [23]. Claims data covering 60 months between May 2013 to April 2018 included 1,163,968 males (age (SD) = 42.5 (13.2) years old) and 1,349,901 females (42.2 (13.1) years old) who were 20 years old or older (but younger than 75 years old). As detailed in Additional File 1 (Appendix 5), we examined 99,821 new users of celecoxib, who used celecoxib after at least 180 days of non-use. The occurrence of peripheral edema was defined by new use of furosemide after at least 180 days of non-use and the index date was the day when the outcome occurred. Daily exposure during an 84-day study period was determined using a 7-day grace period. We selected 311 cases who had both exposed and unexposed days during the study period. The 84-day study period was divided into (M + 1) periods where M = 1, 2, 5, 11, 27 and 83, resulting in periods of 42, 28, 14, 7, 3 and 1 days, respectively. Four different definitions were used to determine exposure in the case and control periods as shown in Table 1. Cases who started celecoxib on the index date were excluded from the analysis since furosemide could have been prescribed for prevention rather than treatment of edema.

The data was analyzed by standard conditional logistic regression (Eq. (2)) and the Mantel-Haenszel method with their 95% CIs, and the weighting method for a binary exposure (Eq. (11)) with 2.5 to 97.5 percentiles estimated by the bootstrap method. We also extended the study period to 168 and 336 days. All the analyses were performed using SAS 9.4. Data and SAS codes to analyze the data are available upon request.

The data for the Australian study on the association between hip fracture and psychoactive medicines [24, 25] were obtained from the Australian Department of Veterans’ Affairs administrative claims database. The study was approved by Department of Defense and Veterans’ Affairs Human Research Ethics (E016–007) and University of South Australia Human Research Ethics (P203/04) where obtaining the informed consent from study subjects was waived for the current study. The cases were 8828 patients aged over 65 years who were hospitalized for hip fracture between 2009 and 2012. The index date for each case was the date of hospitalization. We have re-analyzed the data from the previous case-control study [25] as a case-crossover study with SSRIs as the exposure using the same methods as the Japanese study.

Both of the studies in Japan and Australia were carried out in accordance with the Declaration of Helsinki, and all methods were carried out in accordance with relevant guidelines and regulations in Japan and Australia.

Results

Simulation studies: comparison of methods with or without time-varing confounding

Table 2 and Fig. 3a show OR_SCL, OR_VF, and OR_MH with their 95% CIs, and OR_G with its 2.5 to 97.5 percentiles estimated for the scenario in Fig. 2a for M control periods (1 period = 1 day). The difference between OR_SCL and true RR (4.0) was more than 10% of the true value when M = 6 or > 7 and increased when M increased. When control periods were extended to 10 and 20 cycles (69 and 139 control periods), the estimate (95% CI) of OR_SCL was 7.92 (5.22–12.03) and 8.71 (5.59–13.57), respectively (not shown in Table 2). On the other hand, odds ratios from the remaining 3 methods (OR_VF, OR_MH, and OR_G) were unbiased irrespective of the value of M. The estimate of OR_VF cannot be estimated when M = 7, 10, 14, 17 and 21 because the positivity condition was not satisfied for any data. For example, when M = 7, ∑P(X_i0 = 0, N_exposed = 1) = 0 because X_i0 = 1 in Subgroup H which is only one subgroup where N_exposed = 1 and similarly for N_exposed = 2 or 3, X_i0 was the same for all subgroups.

Table 2 The estimates of OR_SCL, OR_VF, OR_MH, and OR_G from simulated data with a binary exposure only (study period = M + 1 days). Results are shown graphically in Fig. 3a

Fig. 3

Estimates of OR_SCL, OR_VF, OR_MH and OR_G for simulated data in Fig. 1. OR_SCL, OR_VF and OR_MH with their 95% CI, and OR_G with its 2.5–97.5pct for the exposure (exp(β)) for M = 1, 3, 6, 9, 13, 20, 69, and 139, where M is the number of control periods and study period = M + 1 days (1 period = 1 day). (a) shows results for the population in Fig. 2a with one binary exposure variable only. (b) shows results for the population in Fig. 2b with one binary exposure variable and one binary time-varying confounder. OR_SCL: odds ratio by the standard conditional logistic regression (OR_SCL); OR_VF: odds ratio by the Vines and Farrington’s method (OR_VF); OR_MH: odds ratio by the Mantel-Haenszel method (OR_MH); OR_G: odds ratio by the Greenland’s method (OR_G); 95%CI: 95% confidence interval; 2.5–97.5pct: 2.5 to 97.5 percentiles

Table 3 and Appendix Fig. 3b show OR_SCL, OR_VF, and OR_MH with their 95% CIs, and OR_G with its 2.5 to 97.5 percentiles estimated for the time varying confounding scenario in Fig. 2b for M control periods (1 period = 1 day). The difference between OR_SCL for exp(β) and the true RR (4.0) was more than 10% of the true value when M > 5 and increased when M increased. The OR_SCL estimate for exp(β) was larger than the corresponding value in Table 2. When control periods were extended to 10 and 20 cycles, the point estimate of OR_SCL for exp(β) was 9.18 (6.23–13.52) and 10.84 (7.09–16.56), respectively. On the other hand, OR_VF and OR_G for both exp(β) and exp(γ) were in general unbiased, particularly when M > 7 where the estimates of OR_VF and OR_G for exp(β) were within 3% of the true value. The estimate of OR_MH for exp(β) varied between 3.82 (M = 2) and 4.65 (M = 13). Though not shown in Table 3 and Fig. 3bb, OR_MH for exp(β) estimated by ignoring the time-varying confounding was stable but overestimated as 4.67 (17% above the true value). When M = 1, OR_VF was 4.27 for exp(β) (about 7% overestimated) and 1.86 for exp(γ) (7% underestimated). Similarly, when M = 1, OR_G was 4.26 for exp(β) (7% overestimated) and 1.68 for exp(γ) (16% underestimated). When N was increased to 1,000,000, exp(β) and exp(γ) were 3.96 and 2.09 for OR_VF and 3.97 and 2.04 for OR_G when M = 1 (not shown in Table 3).

Table 3 Estimates of OR_SCL, OR_VF, OR_MH, and OR_G from simulated data with a binary exposure and a binary time-varying confounder (study period = M + 1 days). Results are shown graphically in Fig. 3b

Table 4 and Fig. 4a and b show OR_SCL, OR_VF, and OR_MH with their 95% CIs, and OR_G with its 2.5 to 97.5 percentiles for the scenarios in Fig. 2a and b, where the length of the study period was fixed as 22 days but the length of each time period varied between 1 and 11 days. For the scenario in Fig. 2a with one binary exposure variable only, OR_SCL for exp(β) varied between 4.00 and 6.61 when the length of the time period varied, but OR_VF, OR_MH, and OR_G were stable and unbiased (see Table 4 and Fig. 4a). For the scenario in Fig. 2b with one binary exposure variable and one time varying confounder randomly generated in the control periods, the OR_SCL estimate for exp(β) varied between 4.34 (1 period = 11 days) and 6.79 (1 period = 7 days) when M and the length of the time period varied, while OR_VF and OR_G for exp(β) were stable and close to the true value, though the odds ratio was a little overestimated as 4.34 when the time period was 11 days (see Table 4 and Fig. 4b). On the other hand, OR_MH for exp(β) varied 3.54 (1 period = 4 days) and 5.56 (1 period = 11 days). Though not shown in Table 4 and Fig. 4b, OR_MH for exp(β) estimated by ignoring the time-varying confounding was stable but overestimated as 4.67. When M = 1 (1 period = 11 days), OR_VF was 4.31 (8% overestimated) for exp(β) and 1.94 (3% underestimated) for exp(γ). Similarly, when M = 1, OR_G was 4.34 (9% overestimated) for exp(β) and 1.35 (32% underestimated) for exp(γ). When N was increased to 1,000,000, exp(β) and exp(γ) were 3.97 and 1.95 for both OR_VF and OR_G when M = 1 (not shown in Table 4).

Table 4 Estimates of OR_SCL, OR_VF and OR_MH (95% CI), and OR_MH (2.5–97.5pct): Simulated data (study period = 22 days). Results are shown graphically in Fig. 4

Fig. 4

Estimates of OR_SCL, OR_VF, OR_MH and OR_G for simulated data with a fixed study period. OR_SCL, OR_VF and OR_MH with their 95% CI, and OR_G with its 2.5–97.5pct for the exposure (exp(β)) where study period is fixed as 22 days but the number of days in 1 period is varied. When 1 period = 1 or 7 days, OR_VF is not obtainable because the positivity condition is not satisfied. (a) one binary exposure variable only. (b) one binary exposure variable and one binary time-varying confounder. OR_SCL: odds ratio by the standard conditional logistic regression (OR_SCL); OR_VF: odds ratio by the Vines and Farrington’s method (OR_VF); OR_MH: odds ratio by the Mantel-Haenszel method (OR_MH); OR_G: odds ratio by the Greenland’s method (OR_G)

Analyses of case-crossover studies of real-world data

Tables 5 and 6 show the estimates using Exposure Definition I in Table 1 in the Japanese study. In general, OR_SCL was different from OR_G and OR_MH except when M = 1, and OR_SCL estimated from Eq. (2) increased when study period increased: OR_SCL was between 1.91 and 2.82 when study period = 84 days (Table 5), between 2.13 and 4.19 when study period = 168 days (Table 6), and between 3.73 and 7.16 when study period = 336 days (Table 6). The point estimate of OR_G from Eq. (4) was always the same as that of OR_MH, as expected.

Table 5 Estimates of OR_SCL and OR_MH (95% CI), OR_G (2.5–97.5pct): Japanese data on celecoxib-peripheral edema with study period = 84 days and Exposure Definition I

Table 6 Estimates of OR_SCL and OR_MH (95% CI), and OR_G (2.5–97.5pct): Japanese data on celecoxib-peripheral edema with study period = 168 and 336 days and Exposure Definition I

In Appendix Tables 5a and 5b in Additional File 1, the results corresponding to Table 5 and Table 6 for Exposure Definitions II, III and IV (Table 1) are presented. Those results indicated, as in Tables 5 and 6, that OR_SCL varied when M varied as well as when the study period varied, while OR_MH and OR_G were relatively stable. Detailed description for Definition II, III and IV is given in Appendix 5 in Additional File 1.

In the Australian study, after excluding the concordant cases, there were 1316 discordant cases with daily exposed and unexposed periods to SSRIs in the 180 days before the index date. Exposure on the index date was excluded since hip fracture may have occurred the day before admission to hospital.

We used periods of 1, 5, 20, 30, 60 and 90 days with M = 179, 35, 8, 5, 2 and 1, respectively. Exposure within each period was defined using Definition I (Table 1). Estimates of the OR_SCL were biased upwards (Table 7). When M = 1 (exposure period = 90 days), estimates for OR_SCL was identical to OR_MH and OR_G as expected. As in the Japanese study, OR_SCL from Eq. (2) increased with M for M > 1.

Table 7 Estimates of OR_SCL and OR_MH (95% CI), and OR_G (2.5–97.5pct): Australian data on SSRI-hip fracture (Study period = 180 days, Exposure Definition I)

Discussion

Methods to remove BIAS from within-subject exposure dependency with or without time-varying confonders

Using simulated data, we showed that OR_SCL can be biased when there is within-subject exposure dependency but no exposure time trend, except when only one control period is used. When only one control period is used (M = 1), pairwise exchangeability is equivalent to global exchangeability and bias due to within-subject exposure dependency in standard conditional logistic regression does not occur (although bias due to time trends may occur). In Tables 2 and 3, OR_SCL increased with the increase of study period. This observation is similar to that in the previous case-crossover studies, where the odds ratio increased when a longer study period was employed [27,28,29]. In Table 2, OR_VF, OR_MH, and OR_G were unbiased. Of those 3 unbiased estimates, OR_VF may be difficult to calculate when analyzing real-world data because it requires population data (or samples from the population) to estimate the exposure probabilities, unlike OR_MH and OR_G. In addition, probabilities for many exposure permutations must be reliably estimated which may be intractable. For example, in Fig. 1 with 21 control periods, only 8 subgroups A to H were assumed to exist, but in real-world data, many more exposure patterns would occur. It is possible that the positivity condition is not satisfied for certain permutations, and these do not contribute to the estimation of OR_VF. These limitations make OR_VF difficult to use in practical applications.

As Vines and Farrington showed, when no exposure time trend exists and there is one binary exposure variable which is pairwise exchangeable, OR_MH and OR_G are unbiased even when OR_SCL is biased. However, when the model involves a time-varying confounder in addition to the exposure variable, OR_MH may be biased when periods of each subject are further stratified by the time-varying confounder as in Tables 3 and 4. In contrast, the time-varying confounder can be handled by OR_G. Results shown in Tables 3 and 4 indicate that our proposed method is unbiased when both within-subject exposure dependency and time-varying confounding occur at the same time, provided that the time-varying confounder is independent between periods.

In Tables 3 and 4, when M = 1 (i.e., only 1 control period is used), OR_G and OR_VF for exp(β) were overestimated and those for exp(γ) were underestimated, but the estimates approached to the true values when N was increased, suggesting that the deviation from the true values observed when M = 1 was due to random error. Conversely, it is likely that employing a larger number of control periods can produce more precise point estimates of exp(β) and exp(γ).

Real-world data analysis

In the Japanese and Australian studies, we found a discrepancy between OR_SCL and OR_G when more than one control period was used. In the Japanese study, OR_SCL was more than 2 times larger than OR_G when the study period increased. We believe that this discrepancy occurred mainly due to within-subject exposure dependency, which increases as the study period increases. The estimates of OR_G were the same as OR_MH as expected.

In the Australian study, the discrepancy between OR_SCL and OR_G was modest compared with the Japanese study. This was compatible with the finding by Vines and Farrington that the bias due to within-subject exposure dependency is minimal when exp(β) ≈ 1 [15].

Limitations of the current study and future direction

In the current study, our focus was on bias from within-subject exposure dependency. However, biases can occur from other sources as well. First, we have not applied the weighting method when a time trend in the exposure exists [14]. Though it was shown in the methods section that if there is no time trend, probabilities that the individual is exposed (π₁) and unexposed (π₀) at the case period can be estimated without data from population or time-controls, to check whether assumptions of pairwise exchangeability (no exposure time trend) are satisfied, we need data of the population or time-controls. When the proportion of periods exposed in the population or time-controls is roughly constant over study periods, this will support pairwise exchangeability assumptions. When time trend exists, case-time-control design [14] may be used, but more study is needed when both exposure time trend and within-subject exposure dependency exist. Second, a ‘washout period’ has been used in some case-crossover studies [16, 17, 19] to allow for the uncertainty in the optimal length of one period and to reduce within subject exposure dependency. In the current study, when the length of the time period varied, OR_G (and OR_MH) was stable in both the simulation study and real-world data. Since exposure was defined by the last day of the period, any period with 2 or more days was equivalent to using a ‘washout period’ at the beginning of the period. However, much more analyses are needed to examine the need for a ‘washout period’ and to determine the optimal length of one period particularly when within-subject exposure dependency exists. Another bias may occur when the event rate is not constant. For example, the event rate may be particularly high soon after the exposure is started, compared to later after treatment was started. One solution for this problem is to divide the exposed periods into high-risk and low-risk periods and considering this as different levels of exposure. When within-subject exposure dependency exists, this will require weighting for at least 3 exposure levels and the weighting method for binary exposures described in this study should be expanded for multi-level exposures.

Limitations of our study include that we assumed a time-varying confounder with no within-subject dependency. If within-subject dependency exists for a time-varying confounder, the weighting method may need to allow for both the exposure and time-varying confounders. Furthermore, when unmeasured time-varying confounders exist, the results still can be biased even when using our weighting method.

Finally, as mentioned earlier, the variance of OR_G in Eq. (14) is estimated assuming the weight in Eq. (10) is a constant. To allow for uncertainty in estimating the weights in Eq. (10), bootstrapping was used to estimate to 95% CI from 2.5 to 97.5 percentiles of OR_G. However, an analytical method to estimate the variance of OR_G which accounts for the variance of the weights may be useful and be developed in a future study.

Conclusion

Despite autocorrelated exposures being common in pharmacoepidemiology, bias due to within-subject exposure dependency in the case-crossover study has had little attention in the pharmacoepidemiology literature and standard conditional logistic regression has been widely used without assessing whether this bias exists. Although using only one control period can avoid bias due to within-subject exposure dependency, this will reduce the accuracy of the estimate due to random error, as seen in our simulation study (the odds ratio when M = 1 in Tables 3 and 4). To assess for the possibility of bias, we recommend comparing the Mantel-Haenszel odds ratio with the standard conditional logistic regression odds ratio before starting analysis of case-crossover data. If time-varying confounders exist in data set, they may be ignored when comparing two odds ratios to assess the possibility of bias due to within-subject exposure dependency. If a substantial discrepancy is found (e.g., more than a pre-specified threshold such as 5 to 10% of the Mantel-Haenszel odds ratio), standard conditional regression should not be used. Either the Mantel-Haenszel method or our weighting method should be used instead. The weighting method has less bias than the Mantel-Haenszel method when a time-varying confounder exists and in such a case we recommend using our proposed method to estimate OR_G in the analysis of case-crossover data because it removes bias from within-subject exposure dependency and can account for time-varying confounders.

Future research will extend our weighting method to allow for time trends in the exposure, ‘washout periods’, multi-exposure levels which are potentially important when the event rate changes during exposed periods, and analytical method to have the variance of OR_G by making allowance for the uncertainty in estimating the weight.