The Case-time-control Method for Nonbinary Exposures

Abstract

A popular way to reduce confounding in observational studies is to use each study participant as his or her own control. This is possible when both the exposure and the outcome are time varying and have been measured at several time points for each individual. The case-time-control method is a special case, which, under certain assumptions, allows the analyst to control for confounding by time-varying covariates, while controlling for all time-stationary characteristics of the study participants. There are two formulations of the case-time-control method. One formulation requires that the exposure be binary, and the other requires that there be no more than two time points per individual. In this article the author proposes a generalization of the case-time-control method for nonbinary exposures and an arbitrary number of time points. The author derives the asymptotic properties of the resulting estimator and assesses its finite sample properties in a simulation study.

Keywords

case crossover fixed effects regression conditional logistic regression detrending generalized linear model

1. Introduction

A common aim of sociological and epidemiological research is to estimate the causal effect of an exposure (or treatment) on a binary outcome. In observational (nonrandomized) studies, confounding bias is always a concern. If both the exposure and the outcome are time varying and have been measured at several time points for each individual, then confounding by all time-stationary covariates can be eliminated by using each individual as his or her own control. In epidemiological jargon, this strategy is broadly referred to as the case-crossover method (Maclure 1991). A common analytic tool is the conditional logistic regression model, which conditions on the individual and thereby implicitly controls for all individual-specific characteristics that are constant over time (Allison and Christakis 2006).

In many situations the binary outcome is “absorbing”; that is, after it has occurred, the individual is no longer part of the study. The most obvious example is when the outcome is an indicator of death (i.e., 0 for alive, 1 for dead). Often, there may also be time-varying covariates, such as age, income, and socioeconomic status, that we wish to control for in the conditional logistic regression model. However, when at least one of these confounders, such as age, is monotonically increasing with time, it is no longer possible to fit the model. This is because such a confounder perfectly predicts the outcome within the individual, making it impossible to get a within-individual estimate of the coefficient for that covariate (Allison and Christakis 2006).

The case-time-control method (Suissa 1995) has been proposed as a solution to this problem. Briefly, the idea is to let the exposure and outcome switch places in the conditional logistic regression model. Provided that the exposure is not monotonically increasing with time, the model will no longer have convergence problems.

Although both the case-crossover method and the case-time-control method were originally developed in the epidemiological field, they have been recognized as important tools in the sociological field as well (Allison and Christakis 2006), and they have been used in various applications within the sociological field (e.g., Campbell and Lee 2009; Jin and Christakis 2009; Kacperczyk 2013; Sørensen 2007; Winter et al. 2012).

The original formulation of the case-time-control method by Suissa (1995) does not require the exposure to be binary. Indeed, in his application of the method, Suissa also considered dichotomous, trichotomous, and continuous coding of the exposure. However, this formulation works only when the exposure and outcome are observed at exactly two time points for all individuals. Allison and Christakis (2006) used a slightly different formulation of the case-time-control method, which allows an arbitrary number of time points but requires that the exposure be binary. Fafchamps, van der Leij, and Goyal (2010) proposed a related method to control for confounding by time, based on detrending of the exposure. This method works for nonbinary exposures and an arbitrary number of time points, but it is somewhat ad hoc and has not been justified by theoretical arguments.

The aim of this article is to propose a theoretically justified generalization of the case-time-control method for nonbinary exposures and an arbitrary number of time points. The article is structured as follows. Section 2 presents a motivating example with real data, then Section 3 introduces basic notation and definitions and reviews the formulations of the case-time-control method by Suissa (1995) and Allison and Christakis (2006). Section 4 reviews the detrending method proposed by Fafchamps et al. (2010) and shows by example that this method can give considerable bias. Section 5 presents my generalization of the case-time-control method for nonbinary outcomes and an arbitrary number of time points, and Section 6 carries out a small simulation study to assess the finite sample properties of the proposed method. Section 7 returns to the motivating example and reanalyzes this with my proposed generalization of the case-time-control method. Like all statistical methods, the case-time-control method relies on certain assumptions, and Section 8 concludes the article by providing a detailed discussion of these.

2. Motivating Example

As a motivating example I use the data set teenpov borrowed from Allison (2009) and containing information on 1,151 teenage girls who were interviewed annually for five years beginning in 1979. The data set is structured so that one row corresponds to one interview, which makes a total of $1, 151 \times 5 = 5, 755$ rows. Such data are often referred to as “panel data.” I use the following variables: ID (a unique subject-identifier), nonpov (1 if the girl is currently not in poverty according to U.S. federal standards, 0 otherwise), hours (the number of hours currently worked per week), inschool (1 if the girl is currently enrolled in school, 0 otherwise), spouse (1 if the girl is currently living with a spouse, 0 otherwise), age (the girl’s current age), and mother (1 if the girl currently has at least one child, 0 otherwise).

The aim of this article is to investigate how much each additional working hour increases the chance of making a transition from poverty to nonpoverty. With this research question in mind, I restrict attention to those girls who were in poverty at the first interview, and I follow them until the first interview when they were no longer in poverty or until the fifth interview, whichever came first. After this restriction, the data set contains 1,342 interviews or rows from 401 girls. My outcome is nonpov, and my “exposure” is hours. Arguably, the covariates inschool, spouse, age, and mother may be important confounders, which I thus wish to control for in the analysis. I note that all these covariates are time varying; that is, they may change during follow-up within each individual.

As an initial analysis I consider the logistic regression model

\begin{matrix} logit {p (nonpov = 1 | hours, inschool, spouse, age, mother)} \\ = α + β hours + γ_{1} inschool + γ_{2} spouse + γ_{3} age + γ_{4} mother . \end{matrix}

Here, the parameter $β$ is the conditional log-odds ratio between hours and nonpov, given inschool, spouse, age and mother. I fit the logistic regression model with the glm function in R, by typing

> fit <- glm(formula=nonpov~hours+inschool+spouse+ age+mother, family=“binomial”,data=teenpov)

and we obtain the following output:

>summary(fit)

The estimate of $β$ is equal to 0.033, which indicates that each additional working hour increases the odds of making a transition from poverty to nonpoverty, with $\exp (0.033) = 1.03$ units. This association is highly significant (p = $3.06 \times 10^{- 10}$ ), so it is most likely not a chance finding.

However, we may still worry about unmeasured confounding. For instance, it is easy to imagine that unmeasured genetic or lifestyle factors affect both a person’s ability or ambition to work and the person’s ability to make a transition from poverty to nonpoverty. A popular way to deal with unmeasured confounding in sociology and epidemiology is to use the individual as his or her own control; in epidemiological jargon this is referred to as the “case-crossover method.” The method requires repeated measures of the outcome and the exposure, which is exactly what we have in the teenpov data set. Toward this end I consider the conditional logistic regression model

logit {p (nonpov = 1 | i, hours)} = a_{i} + β hours,

where i stands for “individual.” By conditioning on the individual, the model automatically controls for all time-stationary characteristics of the individual, such as genetic make-up (Allison 2009). These covariates are absorbed by the individual-specific intercept $a_{i}$ . Conditional logistic regression is a special case of the broad class of models often referred to as “fixed effects regression models” (Allison 2009). I fit the conditional logistic regression model with the clogit function from the survival package by typing

> library(survival)

> fit <- clogit(formula=nonpov~hours+strata(ID), data=teenpov)

and obtain the following output:

> summary(fit)

coef exp(coef) se(coef) z Pr(>|z|)

hours 0.063238 1.065280 0.009124 6.931 4.17e-12 ***

The conditional logistic regression estimate of $β$ is twice as large as that from the ordinary logistic regression and still highly significant.

However, even though the conditional logistic regression model automatically controls for all time-stationary covariates, it does not automatically control for any time-varying covariates, Thus, we may wish to repeat the conditional logistic regression analysis and control for the variables inschool, spouse, age, and mother, by explicitly adding these as regressors in the model. Unfortunately, though, this is not possible. The reason is that both age and mother are monotonically increasing in time, and the outcome (nonpoverty) is absorbing; that is, it can occur only at the last observed time point for each individual. Thus, when running an analysis that conditions on the individual, like conditional logistic regression, age and mother perfectly predict the outcome, and the model does not converge (Allison and Christakis 2006). This is precisely the problem the case-time-control method intends to solve.

In Section 3 I review the case-time-control method, as formulated by Suissa (1995) and Allison and Christakis (2006). Neither of these formulations works for my data, because my exposure (hours) is nonbinary, and there are more than two time points per individual. In Section 5 I propose a generalization of the case-time-control method that can handle nonbinary exposures with more than two time points, and in Section 7 I reanalyze the data with this generalized method.

3. The Original Formulations of the Case-time-control Method

Let x and y be the time-varying exposure and outcome of interest, respectively. Assume that repeated measures on both x and y are available, and let $x_{ij}$ and $y_{ij}$ be the jth measure of the exposure and outcome, respectively, for individual i, where $i = 1, 2, \dots, n$ and $j = 1, 2, \dots n_{i}$ . Let $t_{ij}$ be the time point at which $x_{ij}$ and $y_{ij}$ are measured for individual i. Note that the time points $t_{i 1}, \dots, t_{i n_{i}}$ do not have to be regularly spaced or the same for all individuals. Assume that the outcome is binary (0/1) and absorbing, so that it can occur only at the last time point; that is, $y_{ij} = 0$ for $j = 1, \dots, n_{i} - 1$ . Let $v_{ij}$ be a vector of time-varying covariates we wish to control for in the analysis, which may contain elements that increase monotonically with time as well as nonincreasing elements. In particular, $v_{ij}$ may include time $t_{ij}$ itself. Finally, I define a time-stationary variable, $g_{i}$ , which is an indicator of whether individual i experiences the outcome at the last time point; that is, $g_{i} = 1$ if $y_{i n_{i}} = 1$ . For all variables, I suppress the indices i and j when not needed. Table 1 illustrates example data for two individuals. The first individual ( $i = 1$ ) is observed at time points 0 and 1 and never experiences the outcome, so that $g = 0$ at $t = 0, 1$ for this individual. The second ( $i = 2$ ) is observed at time points 0, 2, and 4 and experiences the outcome at $t = 4$ , so that $g = 1$ at $t = 0, 2, 4$ for this individual.

Table 1.

Example Data for Two Individuals.

i	j	t	x	y	g
1	1	0	2	0	0
1	2	1	5	0	0
2	1	0	3	0	1
2	2	2	2	0	1
2	3	4	9	1	1

Let $x_{0}$ be an arbitrary exposure level. The target parameter is the log-odds ratio

\log {\frac{p (y = 1 | i, x_{0} + 1, v) p (y = 0 | i, x_{0}, v)}{p (y = 0 | i, x_{0} + 1, v) p (y = 1 | i, x_{0}, v)}},

which measures the association of the outcome with a one-unit increase in the exposure, at given levels of i and v. In principle, this log-odds ratio could depend on both i and v; in statistical jargon, we would then say that there is effect modification by, or interaction with, i and v. I assume, though, that $β$ is constant across levels of i and v.

Ideally we would like to give $β$ a causal interpretation. Toward this end, let $y_{x_{0}}$ be the potential outcome for a given individual when hypothetically assigned a fixed exposure level $x_{0}$ (Pearl 2009; Rubin 1974). In the potential outcome framework we say that we have conditional exchangeability if the potential outcome $y_{x_{0}}$ is conditionally independent of x, given i and v. Conditional exchangeability is a technical condition that holds when i and v are sufficient for confounding control (Pearl 2009). Under conditional exchangeability, $β$ can be interpreted as the causal log-odds ratio

\log {\frac{p (y_{x_{0} + 1} = 1 | i, v) p (y_{x_{0}} = 0 | i, v)}{p (y_{x_{0} + 1} = 0 | i, v) p (y_{x_{0}} = 1 | i, v)}},

which measures the causal effect on the outcome of increasing the exposure with one unit, at a given level of i and v (Pearl 2009).

I consider the conditional (fixed effects) logistic regression model

logit {p (y = 1 | i, x, v)} = a_{i} + β x + γ v .

In this model, all time-stationary covariates are absorbed into the individual-specific intercept $a_{i}$ . For notational convenience I model the effect of v linearly in equation (2), but I emphasize that the linear effect of v can be replaced with a nonlinear function (e.g., splines) in model 2, as well as in all subsequent models.

As mentioned above, standard conditional logistic regression software cannot fit the model in equation (2). This is because time perfectly predicts the outcome within the individual, making it impossible for the estimate of $γ$ to converge (Allison and Christakis 2006). If the exposure is binary (0 or 1), this problem can be solved by assuming a conditional logistic regression model for the exposure, of the form

logit {p (x = 1 | i, y, v)} = b_{i} + β y + δ v .

Model 3 assumes, like model 2, that the log-odds ratio $β$ between the exposure and the outcome is constant across levels of i and v. Because of the symmetry of the odds ratio, it then follows that the parameter $β$ in model 3 is identical to the parameter $β$ in model 2 (e.g., Chen 2007). However, this does not mean that the two models give the same estimate of $β$ . Indeed, whereas model 2 does not converge, model 3 can be fit with conditional logistic regression software without convergence problems, provided that the exposure is not monotonically increasing with time. Let ${\hat{β}}_{x}$ denote the estimate of $β$ on the basis of model 3.

In the formulation of the case-time-control method by Allison and Christakis (2006), $β$ is estimated by fitting model 3. However, Suissa (1995), who originally proposed the method, used a slightly different formulation. Restricting attention to a scenario in which all individuals are observed at exactly two time points, coded as 0 and 1, Suissa proposed using conditional logistic regression software to fit the model

logit {p (t = 1 | i, x, g)} = c_{i} + ψ_{1} x + ψ_{2} xg .

Using Bayes’s rule (together with the definition of the variable g), it can be shown (see Appendix A) that the parameter $ψ_{2}$ in model 4 is equal to the target parameter $β$ , when $v = t$ . This equivalence does not require the exposure to be binary, but it does require that the log-odds ratio between the exposure and the outcome be constant across levels of i and t, as in models 2 and 3. From the equality of $ψ_{2}$ and $β$ , it follows that a consistent estimate of $ψ_{2}$ is also a consistent estimate of $β$ . Let ${\hat{β}}_{t}$ ( $= {\hat{ψ}}_{2}$ ) denote the estimate of $β$ on the basis of model 4.

The two formulations of the case-time-control method, based on models 3 and 4, have different pros and cons. The formulation based on model 3 does not require the time variable to be binary, and it can thus be used when there are more than two time points for some individuals. However, it requires the exposure to be binary, in order for the logistic model for x to make sense. The formulation based on model 4 does not require the exposure to be binary, and it can thus be used for nonbinary exposures as well. However, it requires the time variable to be binary, in order for the logistic model for t to make sense. Also, whereas model 3 can handle additional covariates that increase monotonically with time, model 4 cannot. If both the exposure and the time variable are binary, then it can be shown (see again Appendix A) that models 3 and 4 give the same estimate of $β$ (i.e., ${\hat{β}}_{x} = {\hat{β}}_{t}$ ).

4. Detrending the Exposure

Fafchamps et al. (2010) proposed a related method to control for confounding by time on the basis of a detrending of the exposure. In their proposal, a linear fixed effects regression model is fitted for the exposure, of the form

E (x | i, v) = d_{i} + α v .

The fitted model is used to create residuals of the form $r = x - \hat{α} v$ . Finally, $β$ is estimated by fitting the conditional logistic regression model

logit {p (y = 1 | i, r)} = a_{i} + β r .

Fafchamps et al. (2010) argued informally that this detrending of the exposure should remove any dependence of the exposure on time, so that the resulting estimate of $β$ is not confounded by time.

To gain further insight into this method, suppose we replace the logistic model in equation (2) with the linear fixed effects regression model

E (y | i, x, v) = a_{i} + β x + γ v .

It can be shown (see Appendix B) that the estimate of $β$ obtained from fitting model 7 is identical to the estimate obtained by fitting the fixed effects regression model

E (y | i, r) = a_{i} + β r,

where r is the residual defined above. Thus, detrending the exposure does indeed give the desired result in linear outcome models.

However, there is no theoretical justification for detrending the exposure in logistic outcome models. Indeed, we can easily find examples for which this method gives considerable bias. I give a simulated example here, in which I simulated 1,000 samples of 1,000 individuals each. For each individual, I simulated the normally distributed intercepts $d_{i} ~ N (- 2, 1)$ and $a_{i} ~ N (0, 1)$ . I then simulated repeated measures of a normally distributed exposure and a Bernoulli distributed outcome at five time points $t = 0, 1, 2, 3, 4$ , from models $x | i, t ~ N (d_{i} + 0.4 t, 1)$ and $y | i, t, x ~ Ber {expit (a_{i} + β x + 0.4 t)}$ , respectively, where expit is the inverse of the logit function. The true value of $β$ was equal to 0.5. To mimic an absorbing outcome, I excluded all observations beyond the first time point at which $y = 1$ , for each individual. For each sample I fitted model 5 with $v = t$ and constructed residuals r as described above. Finally, we used these residuals to estimate $β$ in model 6. The mean (over the 1,000 samples) estimate of $β$ was 0.32, which is severely biased. Repeating the simulation with 50 time points $t = 0, 1, \dots, 49$ gave a nearly identical estimate.

5. Generalization of the Case-time-control Method

To generalize the case-time-control method for nonbinary exposures I follow Allison and Christakis (2006) and use a fixed effects regression model for the exposure. This strategy is motivated by the fact that such a model can be fitted for the exposure without convergence problems, provided that the exposure is not monotonically increasing with time (Allison and Christakis 2006). I assume that the exposure has an exponential family distribution

p (x | i, y, v) = \exp {\frac{x θ_{iyv} - A (θ_{iyv})}{ϕ} + k (x, ϕ)}

with canonical parameter $θ_{iyv}$ and dispersion parameter $ϕ$ . I model the conditional mean of x as $η {E (x | i, y, v)} = b_{i} + β^{*} y + δ v$ , where $η$ is either the identity link, the log link, or the logit link. Thus, the model for the exposure becomes a generalized linear model (GLM; McCullagh and Nelder 1989). I assume that $η$ is the canonical link, so that $θ_{iyv} = η {E (x | i, y, v)}$ . Combining this assumption with the fact that $\partial A (θ_{iyv}) / \partial θ_{iyv} = E (x | i, y, v)$ (McCullagh and Nelder 1989) shows that $A (θ_{iyv}) = θ_{iyv}^{2} / 2$ when $η$ is the identity link, $A (θ_{iyv}) = \exp (θ_{iyv})$ when $η$ is the log link, and $A (θ_{iyv}) = \log {1 + \exp (θ_{iyv})}$ when $η$ is the logit link.

As before, let $x_{0}$ be an arbitrary value of x. We see that

\begin{matrix} β = \log {\frac{p (y = 1 | i, x_{0} + 1, v) p (y = 0 | i, x_{0}, v)}{p (y = 0 | i, x_{0} + 1, v) p (y = 1 | i, x_{0}, v)}} \\ = \log {\frac{p (x_{0} + 1 | i, y = 1, v) p (x_{0} | i, y = 0, v)}{p (x_{0} | i, y = 1, v) p (x_{0} + 1 | i, y = 0, v)}} \\ = {\frac{(x_{0} + 1) θ_{i 1 v} - A (θ_{i 1 v})}{ϕ} - k (x_{0} + 1, ϕ)} + {\frac{x_{0} θ_{i 0 v} - A (θ_{i 0 v})}{ϕ} - k (x_{0}, ϕ)} \\ - [{\frac{x_{0} θ_{i 1 v} - A (θ_{i 1 v})}{ϕ} - k (x_{0}, ϕ)} + {\frac{(x_{0} + 1) θ_{i 0 v} - A (θ_{i 0 v})}{ϕ} - k (x_{0} + 1, ϕ)}] \\ = \frac{θ_{i 1 v} - θ_{i 0 v}}{ϕ} \\ = \frac{(b_{i} + β^{*} \times 1 + δ v) - (b_{i} + β^{*} \times 0 + δ v)}{ϕ} \\ = \frac{β^{*}}{ϕ}, \end{matrix}

where the first equality follows from the definition of $β$ , the second equality follows from Bayes’s rule, the third equality follows from the exposure model 9, and the fifth equality follows from the canonical link assumption $θ_{iyv} = η {E (x | i, y, v)}$ together with the assumed mean structure $η {E (x | i, y, v)} = b_{i} + β^{*} y + δ v$ . From equation (10), we observe that the target parameter $β$ is equal to the GLM parameter $β^{*}$ divided by the dispersion parameter $ϕ$ . I thus propose to estimate $β^{*}$ and $ϕ$ , then take the ratio of these to obtain an estimate of $β$ . Let ${\hat{β}}_{GLM}$ denote the resulting estimate of $β$ . If x is binary, $ϕ = 1$ , $k (x, ϕ) = 0$ , and $η$ is the logit link, then the model in equation (9) reduces to the logistic model in equation (3), and ${\hat{β}}_{GLM}$ reduces to ${\hat{β}}_{x}$ . In this sense, my proposal generalizes the “standard” case-time-control method. I proceed by considering nonbinary exposures, for which I assume that $η$ is either the identity link or the log link.

To estimate $(β^{*}, δ)$ , I use conditional generalized estimating equations (CGEE; Goetgeluk and Vansteelandt 2008). This method is applicable when $η$ is the identity link or the log link and provides consistent estimates of $(β^{*}, δ)$ without making any parametric assumptions apart from the mean model $η {E (x | i, y, v)} = b_{i} + β^{*} y + δ v$ . As a by-product, we are able to relax a certain independence assumption, which is needed for consistent estimation when $η$ is the logit link (see Section 8).

To estimate $ϕ$ I proceed as follows. I first define the residual

z (β^{*}, δ) = {\begin{matrix} x - β^{*} y - δ v & if is the identity link \\ x \exp (- β^{*} y - δ v) & if is the log link \end{matrix}

I then estimate $b_{i}$ for each individual as

{\hat{b}}_{i} ({\hat{β}}^{*}, \hat{δ}) = η {\sum_{j = 1}^{n_{i}} z_{ij} ({\hat{β}}^{*}, \hat{δ}) / n_{i}} .

This estimator is motivated by the fact that $E [η^{- 1} {{\hat{b}}_{i} (β^{*}, δ)} | i, y, v] = η^{- 1} (b_{i})$ . I next estimate $θ_{iyv}$ as

{\hat{θ}}_{iyv} ({\hat{β}}^{*}, \hat{δ}) = {\hat{b}}_{i} ({\hat{β}}^{*}, \hat{δ}) + {\hat{β}}^{*} y + \hat{δ} v .

Now, from standard GLM theory we see that $var (x | i, y, v) = ϕ \partial^{2} A (θ_{iyv}) / \partial θ_{iyv}^{2}$ , where $\partial^{2} A (θ_{iyv}) / \partial θ_{iyv}^{2} = 1$ when $η$ is the identity link and $\partial^{2} A (θ_{iyv}) / \partial θ_{iyv}^{2} = \exp (θ_{iyv})$ when $η$ is the log link. Thus, for a given individual i we can construct a moment estimator of $ϕ$ (Pawitan 2001) as

{\hat{ϕ}}_{i} ({\hat{β}}^{*}, \hat{δ}) = \frac{1}{n_{i} - 1} \sum_{j = 1}^{n_{i}} \frac{{[x_{i j} - η^{- 1} {θ_{i y_{i j} v_{i j}} ({\hat{β}}^{*}, \hat{δ})}]}^{2}}{\partial^{2} A {θ_{i y_{i j} v_{i j}} ({\hat{β}}^{*}, \hat{δ})} / \partial θ_{i y_{i j} v_{i j}}^{2} ({\hat{β}}^{*}, \hat{δ})} .

This estimator is not feasible for those individuals with $n_{i} = 1$ . When $η$ is the log link, this estimator is also not feasible for those individuals with $x_{ij} = 0$ for all j, because ${\hat{b}}_{i} ({\hat{β}}^{*}, \hat{δ}) = - \infty$ and $\partial^{2} A {θ_{i y_{ij} v_{ij}} ({\hat{β}}^{*}, \hat{δ})} / \partial θ_{i y_{ij} v_{ij}}^{2} ({\hat{β}}^{*}, \hat{δ}) = 0$ for these individuals. I obtain a “pooled” estimate of $ϕ$ by averaging all feasible individual-specific estimates. Finally, I construct an estimate of $β$ as

{\hat{β}}_{GLM} = \frac{{\hat{β}}^{*}}{\hat{ϕ}} .

To assess the uncertainty in $\hat{β}$ due to sampling variability, it is desirable to derive its asymptotic distribution. I first define $ν = (β^{*}, δ, ϕ, β)$ . The estimator $\hat{ν} = ({\hat{β}}^{*}, \hat{δ}, \hat{ϕ}, {\hat{β}}_{GLM})$ is an M-estimator (Stefanski and Boos 2002) that solves the equation system

\sum_{i = 1}^{n} U_{ν, i} (ν) = \sum_{i = 1}^{n} {\begin{matrix} U_{(β^{*}, δ), i} (β^{*}, δ) \\ U_{ϕ, i} (β^{*}, δ, ϕ) \\ U_{β, i} (ϕ, β) \end{matrix}} = 0,

where $U_{(β^{*}, δ), i} (β^{*}, δ)$ is the CGEE estimating function for $({\hat{β}}^{*}, \hat{δ})$ (Goetgeluk and Vansteelandt 2008), $U_{ϕ, i} (β^{*}, δ, ϕ) = B_{i} {{\hat{ϕ}}_{i} (β^{*}, δ) - ϕ}$ with

B_{i} = {\begin{matrix} 0 & if n_{i} = 1 or (η is the log link and x_{ij} = 0 for all j) \\ 1 & else \end{matrix}},

and $U_{β, i} (ϕ, β) = β^{*} / ϕ - β$ . If $η$ is the identity link and if the underlying assumptions of the case-time-control method hold (see Section 8), then it can be shown that the estimating function $U_{ν, i} (ν)$ is unbiased; that is, $E {U_{ν, i} (ν)} = 0$ (see Appendix C). It follows from standard theory for M-estimators (Stefanski and Boos 2002) that $n^{1 / 2} (\hat{ν} - ν)$ is asymptotically normal with mean 0 and variance given by the “sandwich formula”:

Σ = E' {\frac{\partial U_{ν, i} (ν)}{\partial ν}}^{- 1} var {U_{ν, i} (ν)} E {\frac{\partial U_{ν, i} (ν)}{\partial ν}}^{- 1} .

I obtain a consistent estimate of the variance for $\hat{ν}$ by replacing $ν$ in equation (12) by $\hat{ν}$ and the population moments in equation (12) by their sample counterparts.

If $η$ is the log link, then the estimating function $U_{ϕ, i} (β^{*}, δ, ϕ)$ is biased; that is, $E {U_{ϕ, i} (β^{*}, δ, ϕ)} \neq 0$ . As a consequence, $\hat{ϕ}$ and $\hat{β}$ are biased as well. However, ${\hat{ϕ}}_{i}$ does converge to $ϕ$ as the $n_{i}$ s increase, and thus we would expect that the bias in $\hat{ϕ}$ and $\hat{β}$ goes to 0 as the $n_{i}$ s increase as well. In the next section I confirm this property by simulation.

I end this section with a technical remark. When x has a causal effect on y, the exposure model in equation (9) cannot be interpreted as a data-generating (i.e., causal) model. From a purely mathematical perspective, this is not a problem, because a data-generating model (i.e., a model for the outcome) can always be reformulated (using Bayes’s rule) as a model for the exposure. However, from a subject-matter perspective, we may worry that the postulated exposure model implies a complex or unrealistic (or even impossible) data-generating model. This not the case, though. It follows immediately from results in Chen (2007) that the probability distributions or densities $p (y | i, x = 0, v)$ and $p (x | i, y = 0, v)$ are variation independent of each other, and of the conditional odds ratio between x and y, given i and v. This means that a model for $p (x | i, y, v)$ does not impose any restrictions on $p (y | i, x, v)$ , apart from those implied by the model for the odds ratio. In particular, the exposure model 9 does not rule out (or imply) that the outcome follows a standard logistic model.

6. Simulation

I carried out a small simulation study to assess the finite sample properties of the proposed estimator ${\hat{β}}_{GLM}$ . I first considered a scenario in which $η$ is the identity link. I generated 1,000 samples of 1,000 individuals each. For each individual, I simulated normally distributed intercepts $a_{i} ~ N (- 2, 1)$ and $b_{i} ~ N (0, 1)$ . I then simulated repeated measures of a Bernoulli-distributed outcome and a normally distributed exposure from models $y | i, t ~ Ber {expit (a_{i} + 0.2 t)}$ and $x | i, y, t ~ N (b_{i} + β ϕ y + 0.2 t, ϕ)$ , respectively, where expit is the inverse of the logit function. The true values of $ϕ$ and $β$ were 1.5 and 0.5, respectively. For each individual, the outcome and exposure were generated for $t = 0, 1, \dots, T$ , with T being a common constant for all individuals. To mimic an absorbing outcome, I excluded all observations beyond the first time point at which $y = 1$ , for each individual. In this way the number of time points per individual, $n_{i}$ , becomes a random variable with $\max (n_{i}) = T + 1$ . I first generated data with $T = 1$ . Each sample was analyzed with the method described in Section 5, with $v = t$ . I calculated the mean and standard deviation (over the 1,000 samples) of $\hat{ϕ}$ and ${\hat{β}}_{GLM}$ , together with the mean standard error, as estimated by the “sandwich formula.” We then repeated the procedure for $T = 2$ , $T = 3$ , and $T = 4$ . Table 2 shows the simulation results. We observe that both $\hat{ϕ}$ and ${\hat{β}}_{GLM}$ are virtually unbiased for all values of $\max (n_{i})$ , and their mean standard errors agree well with their standard deviations.

Table 2.

Simulation Results for the Generalized Case-time-control Method When $η$ Is the Identity Link.

	$\hat{ϕ}$ ( $ϕ = 1.5$ )			${\hat{β}}_{GLM}$ ( $β = . 5$ )
$\max (n_{i})$	Mean	SD	Mean SE	Mean	SD	Mean SE
2	1.50	.07	.07	.50	.11	.11
3	1.50	.06	.06	.50	.08	.08
4	1.50	.05	.05	.50	.06	.06
5	1.50	.05	.05	.50	.05	.05

I next considered a scenario in which $η$ is the log link. I repeated the simulation described above but used a different distribution for the exposure. Unless $ϕ = 1$ , the distribution in equation (9) has no simple closed-form expression when $η$ is the log link, which makes simulation difficult. To bypass this problem, I used a common “trick” and instead generated x from a negative binomial distribution with mean equal to $m = \exp (b_{i} + β ϕ y + 0.2 t)$ and variance equal to $m ϕ$ . This distribution has the desired ratio between the variance and the mean, and it converges to the Poisson distribution when $ϕ$ approaches 1 (Pawitan 2001). Table 3 shows the simulation results. Observe that there is a slight bias in both $\hat{ϕ}$ and ${\hat{β}}_{GLM}$ for $max (n_{i}) = 2$ , but this bias goes to 0 as $max (n_{i})$ increases, as expected. All mean standard errors agree well with the corresponding standard deviations.

Table 3.

Simulation Results for the Generalized Case-time-control Method When $η$ Is the Log Link.

	$\hat{ϕ}$ ( $ϕ = 1.5$ )			${\hat{β}}_{GLM}$ ( $β = . 5$ )
max( $n_{i}$ )	Mean	SD	Mean SE	Mean	SD	Mean SE
2	1.39	.06	.06	.54	.08	.08
3	1.44	.05	.05	.52	.05	.05
4	1.47	.05	.05	.51	.04	.04
5	1.50	.04	.04	.50	.03	.03

As a comparison I repeated the simulations above for $T = 1$ but analyzed the data with the conditional logistic regression model 4. I calculated the mean and standard deviation of ${\hat{β}}_{t}$ ( $= {\hat{ψ}}_{2}$ ), together with the mean standard error, as estimated by the inverse Fisher information. Table 4 shows the simulation results. When $η$ is the identity link, ${\hat{β}}_{t}$ in Table 4 is almost unbiased, and its standard deviation is slightly bigger than the standard deviation of ${\hat{β}}_{GLM}$ in the first row in Table 2. When $η$ is the log link, the bias of ${\hat{β}}_{t}$ in Table 4 is similar to the bias of ${\hat{β}}_{GLM}$ in the first row in Table 3; however, the standard deviation of ${\hat{β}}_{t}$ is almost twice the standard deviation of ${\hat{β}}_{GLM}$ .

Table 4.

Simulation Results for the Conditional Logistic Regression Model 4.

	${\hat{β}}_{t}$ ( $β = . 5$ )
$η$	Mean	SD	Mean SE
Identity	.51	.14	.14
Log	.55	.16	.14

Appendix D presents additional simulations in which I used the same setup as described above but with three variations: smaller samples, nonlinear effects of time, and no exposure effect. In these additional simulations we observe similar patterns as in the simulation above, with one notable exception; when there is no exposure effect (i.e., when $β = 0$ ), we observe that the standard deviation of ${\hat{β}}_{t}$ is very similar to the standard deviation of ${\hat{β}}_{GLM}$ . This is an interesting finding for which I currently have no explanation; I recognize this as an interesting topic for future research.

I end this section by following up the technical remark at the end of Section 5. This simulation is based on models for $p (y | i, t)$ and $p (x | i, y, t)$ . When x has a causal effect on y, neither of these models can be interpreted as a data-generating (i.e., causal) model. Thus, we may worry that the postulated models imply a complex or unrealistic data-generating model for the outcome. However, using Bayes’s rule, it is easy to show that the models I have used in the simulation imply the standard logistic regression model $logit {p (y | i, x, t)} = a_{i}^{*} + β x + 0.2 t$ , where $a_{i}^{*} = a_{i} - b_{i} β - β^{2} / (2 ϕ)$ . By marginalizing out y from the postulated exposure model, we see that $p (x | i, t)$ is a mixture of two normal distributions with means equal to $b_{i} + β ϕ + 0.2 t$ and $b_{i} + 0.2 t$ and mixing probabilities equal to $p (y = 1 | i, t) = expit (a_{i} + 0.2 t)$ and $p (y = 0 | i, t) = 1 - expit (a_{i} + 0.2 t)$ , respectively. A simulation from this “data-generating” model is mathematically equivalent to the simulation that we carried out above, but it is less computationally convenient.

7. Motivating Example Revisited

I now return to the motivating example from Section 2. My aim is to estimate the conditional log-odds ratio $β$ between hours and nonpov, given all time-stationary covariates as well as the time-varying covariates inschool, spouse, age, and mother. Following my proposal in Section 5, I used the exposure model 9, with $x =$ hours, $y =$ nonpov, and $v = ($ inschool, spouse, age, mother). I modeled the conditional mean of x as $η {E (x | i, y, v)} = b_{i} + β^{*} y + δ v$ , where I let $η$ be the identity link. I estimated $β^{*}$ and $ϕ$ as described in Section 5, obtaining estimates ${\hat{β}}^{*} = 4.21$ and $\hat{ϕ} = 82.70$ . Finally I obtained the estimate $\hat{β} = {\hat{β}}^{*} / \hat{ϕ} = 0.05$ , which indicates that each additional working hour increases the odds of making a transition from poverty to nonpoverty with $\exp (0.05) = 1.05$ units. This number is close to the estimate I obtained by using conditional logistic regression without controlling for inschool, spouse, age, and mother (see Section 2), which indicates that the confounding by these time-varying covariates may not be very strong in this example. The standard error for $\hat{β}$ was 0.012, which gives a 95 percent confidence interval for $β$ equal to $(0.027, 0.075)$ .

8. Underlying Assumptions of the Case-time-control Method

Like all statistical methods, the case-time-control method relies on certain assumptions. In particular, it relies on the exposure model 9, which assumes that the exposure has an exponential family distribution, with a mean that is linear (after transformation with the link function $η$ ) in i, y, and v. If the model is misspecified, bias may occur. I note that although bias due to model misspecification is a concern, it is not unique for the case-time-control method, but it is a potential problem for all model-based statistical methods. I also note that the model assumptions can be relaxed by, for instance, modeling the effect of v with a polynomial or with splines.

Another issue is that some individuals do not contribute to the parameter estimates when using fixed effects models such as models 3 and 9. To contribute to the estimates, an individual must have variation in x as well as variation in either y or v (or both). When the exposure is binary, it is possible that most individuals are unexposed (i.e., have $x = 0$ ) all through follow-up, in which case most individuals will not contribute to the parameter estimates. This does not violate any assumption of the method per se, but it has important implications for the interpretation of the results. If those few individuals who become exposed at some point (and thus contribute to the estimates) are systematically different (e.g., if they are older) than those individuals who never become exposed, then the results may not be generalizable to the whole population. When the exposure is quantitative, as in my motivating example (Section 2), this problem is less of a concern because most individuals will have variation in the exposure and will thus also contribute to the parameter estimates.

The most important (and subtle) problem with the case-time-control method is that even if the exposure model is correct, certain independence assumptions are required in order for the parameter estimates to be consistent. I begin by defining the vectors $y_{i} = (y_{i 1}, y_{i 2}, \dots, y_{i n_{i}})$ and $v_{i} = (v_{i 1}, v_{i 2}, \dots, v_{i n_{i}})$ . For the logistic exposure model 3, the usual conditional logistic regression likelihood is derived by considering the joint distribution of the exposures $x_{i 1}, x_{i 2}, \dots, x_{i n_{i}}$ , given i, $y_{i}$ , and $v_{i}$ . The derivation proceeds by factorizing this joint distribution into individual-specific contributions, as specified by the individual-level exposure model. However, in order for this factorization to be valid, the following two assumptions are required:

(x_{i 1}, x_{i 2}, \dots, x_{i n_{i}}) are conditionally independent, given (i, y_{i}, v_{i})

and

x_{ij} and (y_{i}, t_{i}) are conditionally independent, given (i, y_{ij}, v_{ij}) .

Assumption (13) states that the exposures at all time points should be conditionally independent, given the individual and the vector of outcomes and the vector of time-varying covariates. Assumption (14) states that the exposure at each time point should be independent of the vector of outcomes and the vector of time-varying covariates, given the individual and the current value of the outcome and the current value of the time-varying covariates. If either of these assumptions is violated, then the usual expression for the conditional logistic regression likelihood is not valid, and the resulting estimates are generally biased (Goetgeluk and Vansteelandt 2008; Zetterqvist et al. 2016). Jensen et al. (2014) carried out a simulation study to investigate the sensitivity of the case-time-control method to assumption (13). They observed some bias when the assumption was violated but noted that the bias was minor for small exposure effects.

Although the reason for requiring assumptions (13) and (14) is rather technical, the assumptions have strong practical implications. Assumption (13) rules out carryover effects from the exposure at one time point to the exposure at later time points, and assumption (14) rules out carryover effects from the outcome at one time point to the exposure at later time points. Such carryover effects are likely to be present in many scenarios. In my motivating example (Section 2), it is clear that both current working hours (exposure) and current poverty status (outcome) may affect working hours at later time points, in which case both assumptions (13) and (14) would be violated to some degree.

However, if we use exposure model 9 with an identity link or log link, only assumption (14) is required. This is because the CGEE technique I use for estimation (see Section 5) does not require that we consider the joint distribution of $x_{i 1}, x_{i 2}, \dots, x_{i n_{i}}$ (Goetgeluk and Vansteelandt 2008; Zetterqvist et al. 2016). Thus, in the reanalysis of the motivating example (Section 7), we may allow current working hours to have an effect on later working hours. However, we must still assume that current poverty status has no effect on the later working hours. It is not obvious how much bias is introduced by violating assumption (14), and I recognize this as an important topic for future research.

In the case-time-control method, conditional logistic regression is used to model the exposure. However, in more standard applications of conditional logistic regression, the outcome is modeled instead. By symmetry, the necessary assumptions are then seen as assumptions (13) and (14), but with the exposure and outcome “swapped,” that is,

(y_{i 1}, y_{i 2}, \dots, y_{i n_{i}}) are conditionally independent, given (i, x_{i}, v_{i})

and

y_{ij} and (x_{i}, v_{i}) are conditionally independent, given (i, x_{ij}, v_{ij}) .

That conditional logistic regression relies on assumptions (15) and (16) is rarely recognized by analysts and rarely stated in standard statistical textbooks. This is unfortunate, because violations of these assumptions may induce bias in the estimates (Sjölander et al. 2012).

9. Discussion

The case-time-control method makes it possible to control for confounding by time while also controlling for all time-stationary characteristics of the study participants. Until now, the method has been restricted to binary exposures or data with only two time points. In this paper I have proposed a generalization of the case-time-control method for nonbinary exposures and an arbitrary number of time points, on the basis of a GLM for the exposure. I have derived the asymptotic properties of the resulting estimator ${\hat{β}}_{GLM}$ , and I have assessed its finite sample properties in a simulation study. The simulation results indicate that the method works well, regardless of the number of time points, when using an identity link in the GLM. When using a log link, ${\hat{β}}_{GLM}$ is somewhat biased for small number of time points. However, for two time points the bias of ${\hat{β}}_{GLM}$ appears to be no larger than the bias of the conventional case-time-control estimator ${\hat{β}}_{t}$ .

Lack of software implementation is often a major obstacle for practitioners who wish to use new methodology. To facilitate the use of the proposed methods, I have written an R function that calculates the estimates of $β^{*}$ , $ϕ$ , and $β$ , together with their “sandwich” standard errors, as described in Section 5. The R function is provided in Appendix E.

Footnotes

Appendix A

Appendix B

Appendix C

Appendix D

Appendix E

Funding

The author disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was funded by the Swedish Research Council (grant 340-2012-6007).

Author Biography

Arvid Sjölander is a statistician and an associate professor in biostatistics in the Department of Medical Epidemiology and Biostatistics at Karolinska Institutet. He received a PhD in causal inference, and he is still working with problems and questions related to causality. His other research topics of interest include attributable fractions, analysis methods for sibling comparison designs, and doubly robust estimation.

References

Allison

Paul D.

2009. Fixed Effects Regression Models: Quantitative Applications in the Social Sciences. Thousand Oaks, CA: Sage.

Allison

Paul D.

Christakis

Nicholas A.

2006. “Fixed-effects Methods for the Analysis of Nonrepeated Events.” Pp. 155–72 in Sociological Methodology, Vol. 36, edited by Stolzenberg

Ross M.

Boston, MA: Blackwell.

Campbell

Cameron D.

Lee

James Z.

2009. “Long-term Mortality Consequences of Childhood Family Context in Liaoning, China, 1749–1909.”Social Science and Medicine 68(9):1641–48.

Chen

Yun H.

2007. “A Semiparametric Odds Ratio Model for Measuring Association.”Biometrics 63(2):413–21.

Fafchamps

Marcel

van der Leij

Marco J.

Goyal

Sanjeev

. 2010. “Matching and Network Effects.”Journal of the European Economic Association 8(1):203–31.

Goetgeluk

Sylvie

Vansteelandt

Stijn

. 2008. “Conditional Generalized Estimating Equations for the Analysis of Clustered and Longitudinal Data.”Biometrics 64(3):772–80.

Jensen

Aksel K. G.

Gerds

Thomas A.

Weeke

Peter

Torp-Pedersen

Christian

Andersen

Per K.

2014. “On the Validity of the Case-time-control Design for Autocorrelated Exposure Histories.”Epidemiology 25(1):110–13.

Jin

Lei

Christakis

Nicholas A.

2009. “Investigating the Mechanism of Marital Mortality Reduction: The Transition to Widowhood and Quality of Health Care.”Demography 46(3):605–25.

Kacperczyk

Aleksandra J.

2013. “Social Influence and Entrepreneurship: The Effect of University Peers on Entrepreneurial Entry.”Organization Science 24(3):664–83.

10.

Maclure

Malcolm

. 1991. “The Case-crossover Design: A Method for Studying Transient Effects on the Risk of Acute Events.”American Journal of Epidemiology 133(2):144–53.

11.

McCullagh

Peter

Nelder

John A.

1989. Generalized Linear Models. 2nd ed. Vol. 37. London: Chapman & Hall.

12.

Pawitan

Yudi

. 2001. In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford, UK: Oxford University Press.

13.

Pearl

Judea

. 2009. Causality: Models, Reasoning, and Inference. 2nd ed. New York: Cambridge University Press.

14.

Rubin

Donald B.

1974. “Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies.”Journal of Educational Psychology 66(5):688–701.

15.

Sjölander

Arvid

Johansson

Anna L. V.

Lundholm

Cecilia

Altman

Daniel

Almqvist

Catarina

Pawitan

Yudi

. 2012. “Analysis of 1:1 Matched Cohort Studies and Twin Studies, with Binary Exposures and Binary Outcomes.”Statistical Science 27(3):395–411.

16.

Sørensen

Jesper B.

2007. “Bureaucracy and Entrepreneurship: Workplace Effects on Entrepreneurial Entry.”Administrative Science Quarterly 52(3):387–412.

17.

Stefanski

Leonard A.

Boos

Dennis D.

2002. “The Calculus of M-estimation.”American Statistician 56(1):29–38.

18.

Suissa

Samy

. 1995. “The Case-time-control Design.”Epidemiology 6(3):248–53.

19.

Winter

Sidney G.

Szulanski

Gabriel

Ringov

Dimo

Jensen

Robert J.

2012. “Reproducing Knowledge: Inaccurate Replication and Failure in Franchise Organizations.”Organization Science 23(3):672–85.

20.

Zetterqvist

Johan

Vansteelandt

Stijn

Pawitan

Yudi

Sjölander

Arvid

. 2016. “Doubly Robust Methods for Handling Confounding by Cluster.”Biostatistics 17(2):264–76.