A General and Robust Estimation Method for the Case-Time-Control Design

Abstract

The case-time-control design is a tool to control for measured, time-varying covariates that increase montonically in time within each subject while also controlling for all unmeasured covariates that are constant within each subject across time. Until recently, the design was restricted to data with only two timepoints and a single binary covariate, or data with a binary exposure. Sjölander (2017) made an important extension that allows for an arbitrary number of timepoints and covariates and a nonbinary exposure. However, his estimation method requires fairly strong model assumptions, and it may create bias if these assumptions are violated. We propose a novel estimation method for the case-time-control design, which to a large extent relaxes the model assumptions in Sjölander. We show in simulations that this estimation method performs well under a range of scenarios and gives consistent estimates when Sjölander’s estimation does not.

Keywords

case-crossover confounding conditional logistic regression conditional maximum likelihood

1. Introduction

Sociological and epidemiological research often aims to estimate the causal effect of an exposure on a binary outcome. One popular way to reduce confounding bias in observational studies is to use each subject as his or her own control, that is, to compare the risk of the outcome at exposed and unexposed timepoints within the same subject. Technically, this is often accomplished by using a conditional logistic regression model. This model conditions on the subject, implicitly controlling for all covariates that are constant within the subject across time, such as genetic makeup (Allison 2009; Allison and Christakis 2006).

Often, scholars will also want to control for measured time-varying covariates in the model, such as age, income, and socioeconomic status. This may be challenging, particularly when (a) some of these covariates are monotonically increasing in time within each subject (e.g., age) and (b) the outcome of interest is “absorbing,” which means the subject leaves the cohort when the outcome occurs (death is the most obvious example). Under these conditions, it it not possible to fit a conditional logistic regression model. This is because a monotonically increasing covariate perfectly predicts an absorbing outcome, so it is impossible to obtain an estimate of the parameter for that covariate (Allison 2009; Allison and Christakis 2006).

The case-time-control design has been proposed as a solution to this problem. Briefly, the design uses a reformulation of the regression model that allows the target parameter to be estimated without convergence problems. There are two common reformulations: one proposed by Suissa (1995) and one proposed by Allison and Christakis (2006). Both of these have important limitations. Suissa’s (1995) reformulation assumes each subject is observed at only two timepoints but allows the exposure variable to be either binary or continuous. Allison and Christakis’s (2006) reformulation allows for more than two timepoints but requires the exposure variable to be binary.

In recent work, Sjölander (2017) proposed an extension of the case-time-control design that allows for an arbitrary number of timepoints, an arbitrary number of (possibly nonbinary) time-varying covariates, and a nonbinary exposure. This extension uses a conditional generalized linear model (GLM) for the exposure. Sjölander (2017) showed that the target parameter (i.e., the exposure effect) can be written as the ratio between the outcome coefficient and the dispersion parameter in this GLM. Sjölander (2017) proposed estimating the outcome coefficient and the dispersion parameter separately using conditional generalized estimating equations (CGEEs) and the method of moments, respectively, and estimating the target parameter as the ratio of these. A major limitation of this two-step estimation method is that it requires the link function in the GLM to be specified. In practice, the analyst often does not know which link function is most appropriate, and an incorrectly specified link function may give biased estimates. Furthermore, the estimation method proposed by Sjölander (2017) only allows for the identity link and the log link; it cannot handle other common link functions in GLMs, such as the inverse or inverse squared.

In this article, we propose a conditional maximum likelihood (CML) estimation method for the case-time-control design that does not suffer from these limitations. This CML method allows for an arbitrary canonical link function in the GLM for the exposure. Furthermore, it does not require the analyst to specify this link function. Thus, this CML method is both more general and more robust than the two-step method proposed by Sjölander (2017).

The article is organized as follows. In Section 2, we introduce notation and assumptions and define the target parameter. In Section 3, we briefly review the two-step method proposed by Sjölander (2017), and in Section 4, we present our CML method. In Section 5, we evaluate the performance of the CML method with a simulation study and compare it to the two-step method proposed by Sjölander (2017). As an illustration of the two-step method, Sjölander (2017) analyzed a real data set containing information on 1,151 teenage girls who were interviewed annually for five years. In Section 6, we reanalyze these data with the CML method and compare the results to those obtained by Sjölander (2017).

2. Notation, Assumptions, and Target Parameter

We assume data have been collected on $n$ subjects at $n_{i}$ time-points for subject $i$ , $i = 1, . . ., n$ . We let $x_{ij}$ and $y_{ij}$ be the exposure and outcome, respectively, for subject $i$ at time $j$ , $j = 1, . . ., n_{i}$ . We assume $y_{ij}$ is binary (0/1) and absorbing (i.e., $y_{ij} = 0$ for $j < n_{i}$ ), but we make no particular assumption about $x_{ij}$ . We let $v_{ij}$ be a column vector of measured covariates for subject $i$ at time $j$ , which we wish to control for in the analysis. Finally, we define $x_{i} = (x_{i 1}, . . ., x_{i n_{i}})$ , $y_{i} = (y_{i 1}, . . ., y_{i n_{i}})$ , and $v_{i} = (v_{i 1}, . . ., v_{i n_{i}})$ .

Most estimation methods for the case-time-control design rely on three independence assumptions (Jensen et al. 2014; Sjölander 2017):

(y_{1}, x_{1}, v_{1}), \dots, (y_{n}, x_{n}, v_{n}) are independent,

(1)

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

(2)

and

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

(3)

Assumption 1 states that the observations from different subjects are independent. This is uncontroversial because it typically holds by study design. Assumptions 2 and 3 are more questionable. Although these are probabilistic assumptions, they have important causal implications. Assumption 2 implies that the exposure at any given timepoint has no direct causal effect on the exposure at future timepoints, and Assumption 3 implies that the exposure at a given timepoint has no direct causal effect on the outcome and covariates at future timepoints, and vice versa. The causal diagram (Pearl 1995, 2009) in Figure 1 illustrates the scenario for two timepoints. In this diagram, an arrow from $x_{i 1}$ to $x_{i 2}$ would violate Assumption 2, and arrows from $x_{i 1}$ into $v_{i 2}$ and $y_{i 2}$ would violate Assumption 3, as would arrows from $v_{i 1}$ and $y_{i 1}$ to $x_{i 2}$ . In real scenarios, these assumptions would often be violated to some extent. However, if the assumptions are approximately valid, then the case-time-control design may give approximately valid inference.

Figure 1.

A causal diagram for which Assumptions 2 and 3 hold.

Note that Sjölander’s (2017) two-step method does not require Assumption 2 to hold. To our knowledge, this is the only proposed estimation method for the case-time-control design that does not require this assumption.

The case-time-control design aims to estimate the odds ratio (OR)

OR (i, x, v) = \frac{p (y_{ij} = 1 | i, x_{ij} = x, v_{ij} = v) p (y_{ij} = 0 | i, x_{ij} = 0, v_{ij} = v)}{p (y_{ij} = 0 | i, x_{ij} = x, v_{ij} = v) p (y_{ij} = 1 | i, x_{ij} = 0, v_{ij} = v)},

which measures the association of the outcome with $x$ units increase in exposure from exposure level $x_{ij} = 0$ for a given subject $i$ at a given covariate level $v_{ij} = v$ . By conditioning on the subject $i$ , the odds ratio implicitly conditions on (i.e., controls for) all covariates that are constant within the subject across time. We assume that the (natural) logarithm of this odds ratio can be parametrized as

\log {OR (i, x, v)} = β^{T} g (v) x,

(4)

where $g (v)$ is some column vector function of $v$ and the target parameter $β$ has the same dimension as $g (v)$ . To allow for a “main effect” of the exposure, one would typically let the first element of $g (v)$ be the constant 1. For instance, if $v_{ij}$ has one element, one could, for instance, define $g (v) = (1, v)^{T}$ and $β = (β_{0}, β_{1})^{T}$ so that $\log {OR (i, x, v)} = β_{0} x + β_{1} xv$ .

A standard way to estimate $β$ is to fit the conditional logistic regression model

logit {p (y_{ij} = 1 | i, x_{ij}, v_{ij})} = a_{i} + β^{T} g (v_{ij}) x_{ij} + γ^{T} h (v_{ij}),

(5)

where $h (v)$ is some column vector function of $v$ and the nuisance parameter $γ$ has the same dimension as $h (v)$ . However, when some components of $v_{ij}$ are monotonically increasing in time, this model will not converge. This is because such covariates perfectly predict the outcome within each subject, thus producing “infinitely large” estimates of $γ$ , which in turn make the obtained estimates of $β$ unreliable (Allison 2009; Allison and Christakis 2006). To solve this problem, one may use the case-time-control design, which reformulates the regression model so that $β$ becomes estimable.

3. The Two-Step Method Proposed by Sjölander (2017)

Sjölander (2017) proposed using a generalized linear model (GLM; McCullagh and Nelder 1989) for the exposure. A GLM assumes that the variable being modeled, in this case the exposure, has a distribution within the exponential dispersion family. This family of distributions is fairly flexible and can accommodate various shapes. It includes, for instance, the Gaussian distribution, the Poisson distribution, the gamma distribution, and the inverse Gaussian distribution. Technically, Sjölander’s (2017) GLM assumes that the conditional distribution for exposure $x_{ij}$ , given subject $i$ , outcome $y_{ij}$ , and covariates $v_{ij}$ , takes the form

p (x_{ij} | i, y_{ij}, v_{ij}) = \exp [\frac{θ_{i} (y_{ij}, v_{ij}) x_{ij} - A {θ_{i} (y_{ij}, v_{ij})}}{ϕ} + k (x_{ij}, ϕ)] .

(6)

In this model, the parameters $θ_{i} (y_{ij}, v_{ij})$ and $ϕ$ are referred to as the canonical parameter and the dispersion parameter, respectively. These generic parameters can be expressed in terms of the standard parameters for each specific member of the exponential dispersion family. For instance, if $x_{ij}$ has a Gaussian distribution with conditional mean $μ_{i} (y_{ij}, v_{ij})$ and constant variance $σ^{2}$ , then $θ_{i} (y_{ij}, v_{ij}) = μ_{i} (y_{ij}, v_{ij})$ and $ϕ = σ^{2}$ , and if $x_{ij}$ has a Poisson distribution with conditional mean $μ_{i} (y_{ij}, v_{ij})$ , then $θ_{i} (y_{ij}, v_{ij}) = \log {μ_{i} (y_{ij}, v_{ij})}$ and $ϕ = 1$ .

The conditional mean of $x_{ij}$ is modeled as $η {μ_{i} (y_{ij}, v_{ij})} = b_{i}^{*} + β^{* T} g (v_{ij}) y_{ij} + δ^{* T} h (v_{ij})$ , where $η$ is an appropriate link function. Sjölander (2017) assumed that $η$ is the canonical link. This means the link function is chosen so that $θ_{i} (y_{ij}, v_{ij}) = η {μ_{i} (y_{ij}, v_{ij})}$ . It follows, for instance, that the canonical link functions for the Gaussian distribution and the Poisson distributions are the identity link and the log link, respectively.

Using Bayes rule, Sjölander (2017) showed that the target parameter $β$ in Model 4 is identical to $β^{*} / ϕ$ in Model 6. This is a key result that allows for estimation of $β$ because the model in Equation 6 can usually be fitted without convergence problems. An important special case occurs when $x_{ij}$ is binary, $ϕ = 1$ , $k (x_{ij}, ϕ) = 0$ , and $η$ is the logit link. In this case, Model 6 reduces to the conditional logistic regression model

logit {p (x_{ij} = 1 | i, y_{ij}, v_{ij})} = b_{i} + β^{T} g (v_{ij}) y_{ij} + δ^{T} h (v_{ij}),

(7)

with $(b_{i}, β, δ) = (b_{i}^{*}, β^{*}, δ^{*})$ . This model can be fitted with standard software provided that $x_{ij}$ is not absorbing. Estimating $β$ by fitting Model 7 corresponds exactly to the case-time-control design as it was formulated by Allison and Christakis (2006). For nonbinary exposures, Sjölander (2017) proposed to first estimate $β^{*}$ and $ϕ$ separately using CGEEs (Goetgeluk and Vansteelandt 2008) and method of moments, respectively, and then estimate $β$ as the ratio of these. He showed that the resulting estimator of $β$ is asymptotically normal, and he derived the asymptotic variance from the “sandwich formula” (Stefanski and Boos 2002).

Notably, CGEEs are restricted to the identity link and log link, and they may give biased estimates if the link function is misspecified. In the next section, we derive an estimation method that allows for an arbitrary canonical link function in the GLM for the exposure and does not require that the analyst specify this link function.

4. The CML Method

Let $x_{i (\cdot)}$ be the order statistic ${x_{i (1)} < . . . < x_{i (n_{i})}}$ . When there are ties in $x_{i}$ , we may rank the tied observations arbitrarily, relative to each other. Under Assumptions 1, 2, and 3 and Model 6, the conditional probability of $x_{i}$ , given ${i, x_{i (\cdot)}, y_{i}, v_{i}}$ , can be written as

p {x_{i} | i, x_{i (\cdot)}, y_{i}, v_{i}} = \frac{\exp [\sum_{j = 1}^{n_{i}} {β^{T} g (v_{ij}) y_{ij} + δ^{T} h (v_{ij})} x_{ij}]}{\sum_{ω \in Ω} \exp [\sum_{j = 1}^{n_{i}} {β^{T} g (v_{ij}) y_{ij} + δ^{T} h (v_{ij})} x_{i ω_{j}}]},

(8)

where $ω = (ω_{1}, . . ., ω_{n_{i}})$ and $Ω$ is the set of all permutations of the integers $(1, . . ., n_{i})$ (see Appendix A). By multiplying over all subjects, we obtain a conditional likelihood $L (β, δ)$ that can be maximized directly to produce an estimate of the target parameter $β$ together with the nuisance parameter $δ$ . It follows from standard likelihood theory that the resulting estimate is asymptotically normal, with a variance-covariance matrix given by the inverse Fisher information matrix. When $x_{ij}$ is binary, $L (β, δ)$ reduces to the usual conditional likelihood for the conditional logistic regression model in Equation 7.

We emphasize that the conditional likelihood $L (β, δ)$ has the same expression regardless of the link function $η$ . This means the analyst does not need to specify the link function; as long as Model 6 holds for some canonical link function $η$ , the CML estimator of $β$ is consistent. This is a major advantage over Sjölander’s (2017) two-step method, which requires the link function to be specified as either the identity or log link and may give biased estimates if it is misspecified.

Kalbfleisch (1978) derived a similar conditional likelihood for nonclustered data in the context of permutation tests. However, in that context, the target parameter is the unscaled parameter $β^{*}$ , whereas in our context, the target parameter is the scaled parameter $β = β^{*} / ϕ$ .

5. Simulation

We carried out a simulation study to assess the finite sample properties of the proposed CML method and compare it with Sjölander’s (2017) two-step method. The code for the simulation is provided in Online Appendix B.

We generated samples of $n$ subjects each, specified in the following. For each subject, we simulated one time-stationary confounder $u_{i}$ from a uniform distribution on the range $(0, 1)$ . For $J$ timepoints, specified in the following, we defined a monotonically increasing confounder as $v_{ij} = j$ , for $j = 1, . . ., J$ . For each subject and timepoint, we simulated a binary outcome from the model $y_{ij} | u_{i}, v_{ij} ~ Ber {expit (- 3 + u_{i} + 0.2 v_{ij})}$ , where expit is the inverse of the logit function. To make the outcome absorbing, we excluded all observations beyond the first timepoint at which $y_{ij} = 1$ for each subject. The number of timepoints per subject, $n_{i}$ , thus becomes a random variable with $\max (n_{i}) = J$ . Finally, for each subject and timepoint, we simulated an exposure from Model 6, with $b_{i}^{*} = u_{i}$ , $β^{*} = β ϕ$ , $β = 0.5$ , and $δ^{*} = 0.2$ . We considered four different exponential family distributions for the exposure: the Gaussian distribution ( $η (μ) = μ$ ), the Poisson distribution ( $η (μ) = \log (μ)$ ), the Gamma distribution ( $η (μ) = - μ^{- 1}$ ), and the inverse Gaussian distribution ( $η (μ) = - μ^{- 2} / 2$ ). For the Poisson distribution, $ϕ$ was set to 1; for the other distributions, $ϕ$ was set to 0.2.

When simulating the exposure from the gamma and inverse Gaussian distributions, we subtracted the constant 2.2 from the confounder $u_{i}$ ; this technicality ensures that the parameters for these distributions are within their range of support. For each exposure distribution, we simulated 10,000 samples. For each sample, we estimated $β$ with three methods: the CML method, the two-step method with an assumed identity link, and the two-step method with an assumed log link. We did not, however, use the two-step method with the log link when the exposure was simulated from the Gaussian distribution because this method requires the exposure to be positive. For each exposure distribution, we calculated the mean and standard deviation (over the 10,000 samples) of the estimates together with the mean standard error as estimated by the “sandwich formula” and inverse Fisher information for the two-step method and the CML method, respectively. The simulation was carried out for the four scenarios: (1) $n = 1, 000$ , $max (n_{i}) = 2$ ; (2) $n = 1, 000$ , $\max (n_{i}) = 5$ ; (3) $n = 300$ , $\max (n_{i}) = 2$ ; and (4) $n = 300$ , $\max (n_{i}) = 5$ .

Table 1 displays the results. When $n = 1, 000$ or $\max (n_{i}) = 5$ (scenarios 1, 2, and 4), the mean standard errors agree well with the corresponding standard deviations for all link functions and estimators. The CML method has a fairly small bias for all link functions, whereas the two-step method has a large bias when the assumed link function is wrong. However, when $n = 300$ and $\max (n_{i}) = 2$ (scenario 3), the CML estimator has a large bias for all link functions except the log link. In this scenario, the mean standard error of the CML estimator appears to underestimate the true standard deviation for the link functions $η (μ) = - μ^{- 1}$ and $η (μ) = - μ^{- 2} / 2$ . For the link function $η (μ) = \log (μ)$ , the mean standard error is much larger than the standard deviation: 18.75 versus 1.65. Notably, the two-step estimator performs well in this scenario when the assumed link function is correct. Note that when the assumed link function is correct, the two-step method has a smaller standard error than the CML method in all scenarios.

Table 1.

Simulation Results When $x_{ij}$ Is Generated from the Endogenous Model 6

		Conditional Maximum Likelihood			Two-Step, Identity Link			Two-Step, Log Link
	$η (μ)$	Mean	SD	Mean SE	Mean	SD	Mean SE	Mean	SD	Mean SE
1	$μ$	.54	.43	.41	.50	.34	.33	—	—	—
	$\log (μ)$	.53	.16	.15	.74	.12	.12	.49	.09	.09
	$- μ^{- 1}$	.55	.50	.49	.82	.48	.47	.76	.51	.50
	$- μ^{- 2} / 2$	.55	.74	.72	.74	.75	.74	.75	.82	.81
2	$μ$	.51	.18	.18	.50	.14	.14	—	—	—
	$\log (μ)$	.51	.06	.06	.71	.05	.04	.48	.03	.03
	$- μ^{- 1}$	.51	.18	.18	.82	.13	.13	.94	.16	.16
	$- μ^{- 2} / 2$	.51	.25	.25	.96	.27	.27	.96	.29	.29
3	$μ$	.70	1.76	1.75	.50	.64	.61	—	—	—
	$\log (μ)$	.77	1.65	18.75	.75	.22	.21	.49	.16	.16
	$- μ^{- 1}$	.69	1.17	.98	.82	.90	.87	.75	.95	.91
	$- μ^{- 2} / 2$	.70	1.56	1.43	.76	1.39	1.34	.78	1.52	1.47
4	$μ$	.53	.33	.33	.51	.26	.26	—	—	—
	$\log (μ)$	.52	.12	.12	.71	.08	.08	.48	.06	.06
	$- μ^{- 1}$	.53	.34	.34	.82	.24	.23	.94	.30	.29
	$- μ^{- 2} / 2$	.52	.46	.46	.95	.50	.49	.94	.53	.52

Note. Target parameter: $β = . 5$ . 1: $n = 1, 000$ , $\max (n_{i}) = 2$ ; 2: $n = 1, 000$ , $\max (n_{i}) = 5$ ; 3: $n = 300$ , $\max (n_{i}) = 2$ ; 4: $n = 300$ , $\max (n_{i}) = 5$ .

In all these simulations, the analysis models were correctly specified: The outcome was generated from an exogenous model, and the exposure was generated from the endogenous Model 6. To assess the performance of the estimators under model misspecification, we modified the simulation scheme so $x_{ij}$ was generated from the exogenous model

p (x_{ij} | i, v_{ij}) = \exp [\frac{θ_{i} (v_{ij}) x_{ij} - A {θ_{i} (v_{ij})}}{ϕ} + k (x_{ij}, ϕ)],

(9)

with $θ_{i} (v_{ij}) = η {μ_{i} (v_{ij})} = b_{i}^{*} + δ^{* T} h (v_{ij})$ , and $y_{ij}$ was generated from the endogenous model $y_{ij} | u_{i}, x_{ij}, v_{ij} ~ Ber {expit (- 3 + β x_{ij} + u_{i} + 0.2 v_{ij})}$ , with the same values of $(b_{i}^{*}, δ^{*}, u_{i}, β)$ as previously described. We considered the same link functions and used the same analysis models as in the aforementioned simulation.

Table 2 displays the results. No scenario or estimator has a substantially larger bias than in Table 1. In some cases, the bias is even smaller (e.g., the CML estimator in scenario 3). These results indicate some robustness of the estimators to model misspecification.

Table 2.

Simulation Results When $x_{ij}$ Is Generated from the Exogenous Model 9

		Conditional Maximum Likelihood			Two-Step, Identity Link			Two-Step, Log Link
	$η (μ)$	Mean	SD	Mean SE	Mean	SD	Mean SE	Mean	SD	Mean SE
1	$μ$	.53	.35	.35	.50	.29	.29	—	—	—
	$\log (μ)$	.51	.11	.11	.67	.09	.09	.45	.08	.07
	$- μ^{- 1}$	.53	.43	.43	.86	.43	.42	.77	.45	.44
	$- μ^{- 2} / 2$	.53	.66	.65	.73	.68	.67	.74	.74	.73
2	$μ$	.50	.16	.16	.50	.13	.13	—	—	—
	$\log (μ)$	.50	.06	.05	.61	.05	.05	.44	.04	.04
	$- μ^{- 1}$	.50	.16	.16	1.14	.13	.13	1.03	.16	.16
	$- μ^{- 2} / 2$	.50	.24	.24	.99	.26	.26	.95	.27	.28
3	$μ$	.60	.72	.67	.50	.54	.53	—	—	—
	$\log (μ)$	.53	.22	.20	.68	.17	.17	.45	.14	.13
	$- μ^{- 1}$	.63	.86	.82	.88	.79	.77	.79	.83	.81
	$- μ^{- 2} / 2$	.65	1.35	1.26	.76	1.26	1.22	.77	1.37	1.34
4	$μ$	.51	.29	.30	.50	.24	.24	—	—	—
	$\log (μ)$	.50	.10	.10	.62	.09	.09	.44	.07	.07
	$- μ^{- 1}$	.52	.31	.30	1.15	.25	.24	1.04	.29	.28
	$- μ^{- 2} / 2$	.52	.44	.44	.99	.48	.47	.95	.50	.50

Note. Target parameter: $β = . 5$ . 1: $n = 1, 000$ , $max (n_{i}) = 2$ ; 2: $n = 1, 000$ , $\max (n_{i}) = 5$ ; 3: $n = 300$ , $\max (n_{i}) = 2$ ; 4: $n = 300$ , $\max (n_{i}) = 5$ .

6. Real Data Analysis

As an illustration, Sjölander (2017) used the data set teenpov, which was borrowed from Allison (2009) and contains information on 1,151 teenage girls who were interviewed annually for five years, beginning in 1979. The data set contains the variables ID (a unique subject-identifier), nonpov. (1 if the girl is currently not in poverty according to U.S. federal standards, 0 else), hours (the number of hours currently worked per week), in school (1 if the girl is currently enrolled in school, 0 else), spouse (1 if the girl is currently living with a spouse, 0 else), age (the girl’s current age), and mother (1 if the girl currently has at least one child, 0 else). To be consistent with previous notation, we rename ID as $i$ ; $i = 1, . . ., 1151$ .

Sjölander (2017) aimed to investigate how much each additional working hour increases the probability of shifting from poverty to nonpoverty. He thus restricted attention to girls who were in poverty at the first interview and followed them until the first interview when they were no longer in poverty or the fifth interview, whichever came first. After this restriction, the data set contains 1,342 interviews from 401 girls, so that $\max (n_{i}) = 5$ .

In a first analysis, Sjölander (2017) fitted the ordinary logistic regression model

\begin{matrix} logit {p (nonpo v_{ij} = 1 | hour s_{ij}, inschoo l_{ij}, spous e_{ij}, ag e_{ij}, mothe r_{ij})} \\ = α^{†} + β^{†} hour s_{ij} + γ_{1}^{†} inschoo l_{ij} + γ_{2}^{†} spous e_{ij} + γ_{3}^{†} ag e_{ij} + γ_{4}^{†} mothe r_{ij} . \end{matrix}

(10)

In Model 10, the parameter $β^{†}$ measures the association between hours and in school, controlling for in school, spouse, age, and mother. The ML estimate of $β^{†}$ is equal to 0.033, which indicates that each additional working hour increases the odds of shifting from poverty to nonpoverty with $\exp (0.033) = 1.03$ units. This causal interpretation may not be justified though because unmeasured confounders could be affecting the association between hours and in school. To reduce confounding bias, one may wish to fit the conditional logistic regression model

\begin{matrix} logit {p (nonpo v_{ij} = 1 | i, hour s_{ij}, inschoo l_{ij}, spous e_{ij}, ag e_{ij}, mothe r_{ij})} \\ = α_{i} + β hour s_{ij} + γ_{1} i nschoo l_{ij} + γ_{2} spous e_{ij} + γ_{3} ag e_{ij} + γ_{4} mothe r_{ij} . \end{matrix}

(11)

By conditioning on the subject-identifier $i$ , Model 11 implicitly conditions on all covariates that are constant within each subject across time. These covariates are, in effect, absorbed into the subject-specific intercept $α_{i}$ . Thus, the parameter $β$ in Model 11 is controlled for all unmeasured time-stationary covariates as well as the measured covariates in school, spouse, age, and mother.

However, Model 11 cannot be fitted to the data because the covariates age and mother increase monotonically in time within each subject and the outcome (nonpoverty) is absorbing. To solve this problem, Sjölander (2017) used the two-step method with an assumed identity link, obtaining an estimate of $β$ equal to 0.051.

To facilitate use of the proposed methods, we wrote an R package, cglm, which is freely available on Cran. We use this package to reconstruct the analysis by Sjölander (2017) and illustrate our novel proposal. We first load the package by typing

> library(cglm)

The package has one single function, cglm, which can be used for both the two-step method and the CML method. To use the two-step method with an assumed identity link and store the results in an object called fit, we type

> fit <- cglm(method=“ts”, formula=hours∼nonpov+inschool+spouse+ age+mother, data=teenpov, id=“ID”, link=“identity”)

The cglm function has five arguments: method, formula, data, id, and link. The method argument specifies whether the two-step method or the CML method should be used: method = “ts” or method = “cml,” respectively. The formula argument specifies the mean structure in the model in the usual R format. The data argument specifies the data frame; we included the teenpov data in the cglm package. The id argument specifies the name of the subject-identifier variable in the data. The link argument specifies the desired link function; this argument is not used when method = “cml”.

We summarize the results by typing

> summary(fit)

Estimate Std. Error z value Pr(>|z|)

nonpov 0.050888 0.010115 5.031 4.88e-07 ***

inschool -0.070630 0.015834 -4.461 8.18e-06 ***

spouse -0.035060 0.037408 -0.937 0.348641

age 0.006511 0.003512 1.854 0.063715.

mother -0.056710 0.015293 -3.708 0.000209 ***

These results are identical to those presented by Sjölander (2017). We caution the reader that the estimated coefficients for in school, spouse, age, and mother in this output, as well as in the following outputs, have no obvious relevance for the research question as these measure the association between the measured covariates and the exposure (scaled by the dispersion parameter).

To use the two-step method with an assumed log link, we type

> fit <- cglm(method=“ts”, formula=hours∼nonpov+inschool+spouse+age+mother, data=teenpov, id=“ID”, link=“log”)

> summary(fit)

Estimate Std. Error z value Pr(>|z|)

nonpov 0.019020 0.010531 1.806 0.07091.

inschool -0.039352 0.013810 -2.850 0.00438 **

spouse -0.027948 0.039615 -0.705 0.48051

age 0.016113 0.006397 2.519 0.01177 *

mother -0.056161 0.027556 -2.038 0.04154 *

The estimate of $β$ is now equal to 0.019. Thus, for these data, the two-step method seems fairly sensitive to the choice of link function.

To use the CML method, we type

> fit <- cglm(method=“cml”, formula=hours∼nonpov+inschool+spouse+age+mother, data=teenpov, id=“ID”)

> summary(fit)

Estimate Std. Error z value Pr(>|z|)

nonpov 0.046174 0.011689 3.950 7.81e-05 ***

inschool -0.042831 0.012195 -3.512 0.000444 ***

spouse -0.037398 0.023299 -1.605 0.108459

age 0.008594 0.004281 2.007 0.044701 *

mother -0.047147 0.019281 -2.445 0.014475 *

The CML estimate of $β$ is equal to 0.046, which is close to the two-step estimate with an assumed identity link. This indicates that each additional working hour increases the odds of shifting from poverty to nonpoverty with $\exp (0.046) = 1.047$ units. The p value is close to 0, which indicates this finding is statistically significant.

In the previous reanalysis, we assumed that the exposure-outcome odds ratio is constant across levels of covariates (e.g., $g (v) = 1$ in Model 4). As an illustration, we fitted a slightly more elaborate model, allowing for an interaction between nonpov. and spouse:

> fit <- cglm(method=“cml”, formula=hours∼nonpov*spouse+inschool+age+mother, data=teenpov, id=“ID”)

> summary(fit)

Estimate Std. Error z value Pr(>|z|)

nonpov 0.048108 0.012315 3.906 9.37e-05 ***

spouse -0.021928 0.033308 -0.658 0.510318

inschool -0.042947 0.012254 -3.505 0.000457 ***

age 0.008501 0.004300 1.977 0.048066 *

mother -0.046529 0.019442 -2.393 0.016700 *

nonpov:spouse -0.023196 0.041239 -0.562 0.573798

The interaction term on the last row is nonsignificant, and the main effect on the first row is almost identical to the main effect in the aforementioned simpler analysis.

7. Discussion

The case-time-control design makes it possible to control for measured, time-varying covariates that increase monotonically in time within each subject while also controlling for all unmeasured covariates that are constant within each subject across time. Until recently, the design was restricted to data with only two timepoints and a single binary covariate, or data with a binary exposure. Sjölander (2017) made an important extension that allows for an arbitrary number of timepoints and covariates and a nonbinary exposure. However, his two-step estimation method requires specification of the link function in the GLM for the exposure and is restricted to the identity link and the log link.

In this article, we proposed a novel CML estimation method for the case-time-control design. This method allows for an arbitrary canonical link function and does not require the analyst to specify the link function. Our simulations show that the CML method works well and delivers a consistent estimate of the target parameter for a wide range of link functions. In contrast, we have shown that Sjölander’s (2017) two-step estimation method may be highly sensitive to the choice of link function and may have a large bias if this is misspecified.

Despite the strong dependency on the link function, we note three potential advantages of the two-step method. First, our simulations show that the two-step estimator is more efficient than the CML estimator when the link function is correct. This is not surprising given that the CML method makes weaker model assumptions than does the two-step method and the CML method conditions on the order statistic $x_{i (\cdot)}$ , which may contain information about $β$ . Second, when samples are small (i.e., few subjects and timepoints), we find a large bias for the CML estimator but not for the two-step estimator. This indicates that the two-step estimator may perform better in small samples than the CML estimator provided the assumed link function is correct.

A third, more subtle advantage of the two-step method is that it does not require the independence Assumption 2 to hold. Whereas the CML estimator requires that the exposure at a given timepoint has no direct causal effect on the exposure at future timepoints, the two-step estimator does not. This is potentially an important advantage as exposure-to-exposure causal effects are likely present in many scenarios. Unfortunately, it is difficult to study violations of this assumption in isolation: Both the two-step method and the CML method require that the exposure follows the GLM in Equation 6 at each timepoint, and it is difficult to see how one could make this model hold while violating the independence Assumption 2. We recognize this as an important topic for future research.

Supplemental Material

SM862259_Supplemental_Appendix_B – Supplemental material for A General and Robust Estimation Method for the Case-Time-Control Design

Supplemental material, SM862259_Supplemental_Appendix_B for A General and Robust Estimation Method for the Case-Time-Control Design by Arvid Sjölander and Yang Ning in Sociological Methodology

Footnotes

Appendix A

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Funding was provided by the Swedish Research Council (grant no. 2016-01267).

Author Biographies

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

Yang Ning is a statistician and assistant professor in the Department of Statistical Science at Cornell University. He received his PhD in biostatistics from the Johns Hopkins University. His research interests include high-dimensional statistics, causal inference, and statistical inference under nonstandard conditions.

References

Allison

Paul D.

2009. Fixed Effects Regression Models, Quantitative Applications in the Social Sciences 160. Los Angeles: SAGE.

Allison

Paul D.

Christakis

Nicholas A.

2006. “Fixed-Effects Methods for the Analysis of Nonrepeated Events.” Sociological Methodology 36:155–72.

Goetgeluk

Sylvie

Vansteelandt

Stijn

. 2008. “Conditional Generalized Estimating Equations for the Analysis of Clustered and Longitudinal Data.” Biometrics 64: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:110–13.

Kalbfleisch

John D.

1978. “Likelihood Methods and Nonparametric Tests.” Journal of the American Statistical Association 73:167–70.

McCullagh

Peter

Nelder

John A.

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

Pearl

Judea

. 1995. “Causal Diagrams for Empirical Research.” Biometrika 82:669–88.

Pearl

Judea

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

Sjölander

Arvid

. 2017. “The Case-Time-Control Method for Nonbinary Exposures.” Sociological Methodology 47:182–211.

10.

Stefanski

Leonard A.

Boos

Dennis D.

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

11.

Suissa

Samy

. 1995. “The Case-Time-Control Design.” Epidemiology 6:248–53.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB