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

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. Supplementary Information The online version contains supplementary material available at 10.1186/s12874-021-01408-5.


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 withinsubject 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 withinsubject 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 realworld 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 exposureoutcome 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 casecrossover 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]: 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 where N exposed is the total number of exposed (case and control) periods, given by N exposed = M m=0 X im and κ 1 P X i0 = 1, N exposed = k is the sum of probabilities for all the permutations of k exposed and (M + 1-k) unexposed periods where X i0 = 1, and κ 0 P X i0 = 0, N exposed = k 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 .
Vines and Farrington showed that the likelihoods in Eqs. (1) and (2) (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 (π 0 ) and exposed (π 1 ) 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 π 0 , π 1 , and the relative exposure probability π 10 = π 1 / π 0 are estimated, the following likelihood for casecrossover studies proposed by Greenland [26] can be used to obtain an unbiased estimate of exp(β), defined as OR G : 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 /π 0 . 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 where x ic = 1 and π c0 = π 10 = π 1 /π 0 when the individual i was exposed at the case period and x ic = 0 and π c0 = π 00 = π 0 /π 0 = 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: In Eq. (5), X 0 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 0 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 Using conditional probabilities, pairwise exchangeability, P{X i0 = 1, X im = 0} = P(X i0 = 0, X im = 1} can be rewritten as: where π 0 = P(X 0 = 0) and π 1 = P(X 0 = 1). When both sides of Eq. (6) are summed up over M control periods (m = 1, 2, −-, M), we obtain: The quantity M m=1 P( X m = 0|X 0 = 1) can be estimated by the average number of unexposed control periods (where X m = 0) in those exposed at the case period (X 0 = 1). Similarly, the quantity M m=1 P( X m = 1|X 0 = 0) can be estimated as the average number of exposed control periods (X m = 1) in those unexposed at the case period (X 0 = 0). This average, defined as PT 10 and PT 01 , respectively, can be written as: 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 1 is the number of discordant exposed cases with at least one unexposed control period and a 0 is the number of discordant unexposed cases with at least one exposed control period. When M m=1 P( X m = 0|X 0 = 1) and M m=1 P( X m = 1|X 0 = 0) in Eq. (7) are substituted by PT 10 and PT 01 , respectively, we obtain: 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: where m 1 i is the number of exposed periods and m 0 i is the number of the unexposed periods (including both case and control periods) in case i and π 10 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: where w ij = w i1 and w ij = w i0 when j-th period is exposed and unexposed, respectively, in case i.
Using a 1 and a 0 , Eq. (4) can be rewritten as Equation (12) gives (see Additional File 1, Appendix 2) the following maximum likelihood estimate for OR G : When π 10 estimated by Eq. (9) is considered as a constant, the variance is given by: 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: 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 (11) 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 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 1 ) and during unexposed periods (f 0 ) 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). 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 0 ) 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 0 when unexposed and N RR r 0 when exposed 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 = 0.2 and f 1 = 0.4 where f 0 and f 1 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 and N RR RR z r 0 f 1 , 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.

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 reanalyzed 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. 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  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 Table 1 Exposure Definitions Last day = exposure status of the period was the exposure status on the last day; Half or more days = exposure status of the period is defined as 'exposed' if at least half of days during the period was exposed and 'unexposed' otherwise; Any 1 day = exposure status of the period is defined as 'exposed' if at least 1 day during the period is exposed and 'unexposed' otherwise  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 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 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). Tables 5 and 6 show the estimates using Exposure Definition I in Table 1 in the Japanese study. In general, OR SCL 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 M is the number of control periods and study period = (M + 1) days. 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.

Analyses of case-crossover studies of real-world data
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.

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 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 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 M is the number of control periods and study period = (M + 1) days.  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 timevarying 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 withinsubject 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 (π 1 ) and unexposed (π 0 ) 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 timevarying 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 timevarying 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 prespecified 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.