Sage Journals: Discover world-class research

Abstract

Models for time-to-event data are based on transition rates between states, and to define such hazards of experiencing an event, the time scale over which the process evolves needs to be identified. In many applications, however, more than one time scale might be of importance. Here, we demonstrate how to model a hazard jointly over two time dimensions. The model assumes a smooth bivariate hazard function, and the function is estimated by two-dimensional $P$ -splines. We provide an R package for the analysis of event history data with two time scales. As an example, we model transitions from cohabitation to marriage or separation simultaneously over the age of the individual and the duration of the cohabitation. We use data from the German Family Panel (pairfam) and demonstrate that considering the two time scales as equally important provides additional insights about the transition from cohabitation to marriage or separation.

Keywords

Time-to-event data bivariate hazard smoothing P-splines cohabitation marriage separation

Introduction

Analysis of time-to-event data focuses on the timing of successive transitions between states of interest and the effect of covariates on the occurrence of these transitions. The life course can be viewed as a sequence of such transitions or events, and a common statistical framework to analyze how life courses evolve is a stochastic process with a continuous time scale and a finite number of states (Blossfeld et al. 1989), that is, a multi-state model. In its most reduced form, the former simplifies to survival analysis. Only one event is considered: the process starts at birth and is ended by the event “death” and the time scale is the age of the individual. When modeling transitions from one state to another in the life course, hazard models are the most common approach. Hazard functions capture the time-dependent instantaneous risk of an event, given the transition has not yet occurred, and by definition depend on the time scale of the process.

Time is, therefore, key in understanding any event history analysis, and each time scale is defined with respect to a reference event, for example, birth (Hobcraft and Murphy 1986). In demographic analyses, the most prominent time scale is age, which is time since birth. This is because many of the life transitions, such as leaving parental home, completing education, entering a first union or birth of the first child, are thought to be age-dependent. Demographers refer to statistical regularities in the timing of demographic events in a population as “age norm” (Settersten 2003). Consequently, age-specific transition rates are the most frequent choice.

But different time scales can be relevant when analyzing the occurrence of an event like, for example, marriage or divorce. It has long been argued that other time dimensions, such as the duration since entry into the current state, are also of interest (Hobcraft and Murphy 1986). For example, if the timing of the first child birth within a marital union is of interest, the hazard of first birth strongly depends on the age of the mother, as women’s fertility declines with increasing age. However, the duration of the marriage is of equal importance, as couples might base the decision of having a child on the time elapsed since their marriage. The duration of the marriage, therefore, represents another time scale over which the hazard of marital child birth changes.

Standard hazard models require the choice of one time scale over which the risk of observing an event evolves (Klein and Moeschberger 2005). This choice is usually based on previous knowledge of the phenomenon under study. While in demography, age-specific rates are commonly preferred, in medical and epidemiological studies, when disease progression is analyzed, time since diagnosis or disease onset might be the first choice. If additional time scales are present and acknowledged, a common strategy is to incorporate the second time dimension as a categorized fixed or time-varying covariate in the model. Alternatively, a more flexible approach often used to deal with non-proportionality of the hazards is to account for time-varying effects of a fixed covariate, for example, by estimating a time-varying effect of the age at entry into the risk set (Bender et al. 2018).

In demography, there have been attempts to model the hazard of an event directly as a function of more than one time dimension. An early attempt to model the hazard of marital childbearing by age and duration of marriage comes from Page (1977), who models fertility rates as a function of age, cohort, and marriage duration. The Page model is rather inflexible in the way the hazard depends on these time scales, and it requires a parametric formulation of the effect of each time scale on the hazard. Other examples of models that consider hazard functions simultaneously over two time axes can be found in the medical statistical literature (Batyrbekova et al. 2022; Bower et al. 2022; Efron 2002; Iacobelli and Carstensen 2013; Rebora et al. 2015), but the issue has rarely been explored in the social sciences.

With this paper, we aim to shed light on the overlooked aspect of multiple time scales in event history analysis within a field of application where it has not received thorough examination. We review some approaches to model hazards over two time scales simultaneously, which differ in how flexibly the effect of the second time scale is modeled. The models make only modest assumptions about the bivariate hazard, namely that the function evolves smoothly over both time domains. We estimate the models by $P$ -splines.

To illustrate the approach, we study transitions from unmarried cohabitation to marriage or to separation by age of the individual and duration of the cohabitation. Age certainly has a major role in determining the transitions from cohabitation to a subsequent state; however, previous research has favored duration of the cohabitation as the main time scale for such transitions (Di Giulio et al. 2019; Hiekel et al. 2015). We argue that the two time dimensions will likely interact with each other and that the hazard of marrying or separating is determined simultaneously by both time scales. We explore these questions using data from the German Family Panel (pairfam) (Brüderl et al. 2020a).

This study contributes to the literature on modeling time-to-event socio-demographic data in several ways. First, we provide theoretical and practical reasoning to consider more than one time scale for the analysis of time-to-event data, with a particular focus on demographic processes. Second, we review the $P$ -splines approach to smooth hazards with two time scales. Finally, we present an application in family demography, and we demonstrate that with our flexible model, we are able to discover previously overlooked patterns in the transition rates out of unmarried cohabiting unions.

In Section “Temporal Change: Time Scales, Covariates, Effects,” we first discuss time scales in event history analysis and how additional time scales can enter as covariates in the model. Then, we describe the two-dimensional $P$ -spline model for the hazard of an event. Section “ Transitions From Cohabitation to Marriage or Separation” is devoted to the application. We first discuss the demographic background and the data, and then show the results of our two-way model. A discussion concludes the paper.

Temporal Change: Time Scales, Covariates, Effects

Social sciences are concerned, among other things, with the study of temporal change. On the individual level, change often is marked by the occurrence of events, such as obtaining an educational degree, moving out of parental home, marriage, birth of a child, divorce, entering retirement, and ultimately death—just to name a few out of a long list. The technical vehicle to study such histories of events is stochastic processes, that is, random variables indexed by time that describe the occurrence of the event of interest probabilistically.

The time scale over which the process is modeled commonly is obvious—the age of the individual or time since marriage—but temporal change often occurs along several time scales simultaneously. Individuals change their attitudes and behavior with age, but also macro environment changes with calendar time, which influence the occurrence of the event of interest. The concept of individual life courses evolving along several time scales is long known, and portraying processes along two time scales is familiar to social scientists ever since Lexis (1875). If more than one time scale is relevant, then all should enter the analysis, and their importance should be assessed properly.

When the occurrence of events is studied, then hazard-based models are the most common choice. For a single time scale $t$ , the hazard $λ (t)$ is the instantaneous risk of the event occurring at time $t$ , given it has not occurred yet

λ (t) = lim_{ν ↓ 0} \frac{P (e v e n t i n [t, t + ν) | n o e v e n t b e f o r e t)}{ν} .

The choice of the time scale in such models can affect how the risk factors are associated with the rate of the event (Wolkewitz et al. 2016). How to optimally select which of the time scales to be used (in case only one is to be considered) in a time-to-event model has been a topic of discussion in the medical statistical literature (Pencina et al. 2007; Thièbaut and Bènichou 2004), and will not be discussed here.

Time scales differ only in their origin, and they move at the same speed, so that an increment of one time unit on one scale corresponds to the same increment on the other scale. Because of the same speed, when portrayed in a so-called Lexis diagram, individuals move along a diagonal line with slope 1 (Keiding 1990) (see Figure 1(a)). In the application presented in Section “Transitions From Cohabitation to Marriage or Separation,” we will study transitions out of unmarried cohabitation (separation or marriage) and the two time scales will be age (denoted by $t$ ) and duration of the cohabitation (denoted by $s$ ). In the following (and also in Figure 1(a)), for the sake of clarity, we will repeatedly refer to the application, but the concepts apply in general.

Figure 1.

Eight individual trajectories for cohabitation until marriage, separation, or censoring. The panel on the left shows the trajectories in the $(t, s)$ -plane. The panel on the right shows how these trajectories are mapped into vertical trajectories in the $(u, s)$ -plane, where the events are scattered in the full positive plane. The orange points indicate transitions to marriage, the light blue triangles indicate transitions to separation, and the black squares indicate right-censored observations. Individuals may be observed from the beginning of cohabitation (A) or they may enter the sample after cohabiting already for some time (B), which implies left truncation.

The two time scales hazard $λ (t, s)$ is defined as

λ (t, s) = lim_{ν ↓ 0} \frac{P (event \in {t + ν, s + ν} | no event before (t, s))}{ν},

(1)

and it is the instantaneous risk of experiencing an event at age $t$ and duration of the cohabitation $s$ , given that no event happened before $(t, s)$ . The set ${t + ν, s + ν}$ denotes the 45-degree line induced by the same increment along $t$ and $s$ .

In the example, time scale $s$ starts later (at the beginning of cohabitation) than time scale $t$ (which starts at birth). This implies that $t > s$ (it is impossible to be cohabiting for longer than being alive) and consequently individual life lines will be found only in the lower sub-triangle of the positive quadrant in the $(t, s)$ -plane. We will adhere to this notational principle (i.e., $t > s$ ) in the paper. The second time scale $s$ often is kicked off by an event itself. If $t$ is calendar time, then birth starts time scale age $s$ . If $t$ is age, then moving in together starts the duration of cohabitation $s$ ; and so forth. (When embedded in a multistate model that also spans the process up to the start of the later time scale, then $s$ is commonly referred to as clock-reset time scale and can be used to check the Markov assumption in the multistate model, see Tassistro et al. 2019.)

The point in time $t$ when the second scale commences ( $s = 0$ ) gives the difference in the origins of the two time scales, and we denote it by $u = t - s$ . Other than $t$ and $s$ , which continuously change along the 45 $\circ$ lines, the value of $u$ is fixed for each individual (but in general differs across individuals). In the first example above, $u$ is the date of birth, in the second example, $u$ is the age at which the individual enters cohabitation. The linear equation $u = t - s$ is the notorious cohort = period $-$ age relationship troubling some of the age period cohort analysis: Knowing two of the quantities determines the third.

In the cohabitation example, employing $(u, s)$ instead of $(t, s)$ makes individuals move along vertical lines in their “cohabitation histories,” see Figure 1(b). In this way, the observations are, in principle, located in the full positive plane $(u, s)$ , which can make some modeling approaches more handy due to the less uncommon domain. Note, however, the transformation from $(t, s)$ to $(u, s)$ affects only the coordinates of the data points, and not the level of the hazard. So we have

\overset{˘}{λ} (u, s) \equiv λ (t, s),

(2)

where $\overset{˘}{λ} (u, s)$ denotes the value of the hazard when expressed in $(u, s)$ coordinates. Using $t = u + s$ , the estimated hazard can be transformed back to the $(t, s)$ domain unambiguously.

The elaboration above demonstrates that for estimating a hazard along two time scales, $λ (t, s)$ , one can pursue modeling in the $(u, s)$ -plane using equivalence (2). In the latter case, $u$ is a fixed covariate, commonly produced as an outcome of an upstream process, and the importance (or not) of the second time scale $t$ will manifest in how far the covariate $u$ modulates a single-time-scale hazard $\overset{˘}{λ} (s)$ .

In the simplest case, when only one time scale matters, then $\overset{˘}{λ} (u, s) = λ_{0} (s)$ . Potential impact of $u$ (and hence of the second time scale $t$ ) can be considered in a hierarchy of complexity, see Table 1. We formulate all models for the log-hazard $\ln \overset{˘}{λ} (u, s)$ , a common strategy in hazard modeling. In all models, we want to avoid parametric specifications and only make the modest assumption that the functions $θ$ , $ϕ$ , and $η$ , respectively, are smooth.

Table 1.

Models for $\ln \overset{˘}{λ} (u, s)$ .

Model hierarchy $\ln \overset{˘}{λ} (u, s) = \overset{˘}{η} (u, s)$
(A)	$\ln \overset{˘}{λ} (u, s) = \overset{˘}{η} (u, s) = θ (s) + ϕ (u)$	log-additive
(B)	$\ln \overset{˘}{λ} (u, s) = \overset{˘}{η} (u, s) = θ (s) + ϕ (s) u$	time-varying coefficient
(C)	$\ln \overset{˘}{λ} (u, s) = \overset{˘}{η} (u, s)$	full 2D interaction

Note. 2D = two-dimensional.

In Model (A), the same baseline shape $θ (s)$ prevails over $u$ , only the level of the hazard is modulated by $ϕ (u)$ . This is a proportional hazards model; the simplest case is a linear $ϕ (u) = β u$ . Should $ϕ (u)$ be constant, the model reduces to a single time scale hazard over $s$ . Model (B), the varying-coefficient model, allows “interaction” between $s$ and $u$ , but in a restricted way. The impact of the covariate $u$ is no longer constant over $s$ , but the effect of $u$ is time-varying (along scale $s$ ). The way the effect changes with $s$ , that is, $ϕ (s)$ is identical across the values of $u$ though. Model (C) releases this assumption and estimates the log-hazard $η (u, s)$ as a smooth surface over the positive $(u, s)$ plane so that $u$ and $s$ can “interact” flexibly.

For estimating Models (A) to (C), we will employ univariate or bivariate $P$ -splines (Eilers and Marx 2021). The technical details will be described in Section “Smoothing Bivariate Hazard rates with P-splines.” In Section “Transitions From Cohabitation to Marriage or Separation,” we will demonstrate the approach for the application. We also graphically illustrate the differences in the three models (A)–(C).

Alternative approaches to model and estimate hazards as functions over two time dimensions have been described in the literature. The two-way model proposed by Efron (2002) expresses the log-hazard as the sum of two baseline hazards, polynomials in this case, each estimated over one of the two time scales. More versatile versions of this model include flexible parametric functions of the two time scales (Batyrbekova et al. 2022; Bower et al. 2022; Iacobelli and Carstensen 2013). A full interaction model, although in principle possible, to our knowledge has not been implemented yet. Log-additive and varying-coefficient specifications have also been proposed in Bender et al. (2018), see also Bender and Scheipl (2018).

Smoothing Bivariate Hazard Rates With P-splines

Smoothing a univariate (single time scale) hazard $λ (t)$ can be conveniently done using $P$ -splines (Eilers and Marx 1996). The time axis is divided into (narrow) bins $I_{k}$ , of equal width, and within each bin, the number of observed events $y_{k}$ and the time at-risk (exposure) $r_{k}$ is determined from the sample. If the hazard in bin $k$ is $λ_{k}$ , then the likelihood for $λ = (λ_{1}, \dots, λ_{K})^{T}$ is identical to a Poisson likelihood when the counts $y_{k}$ are realizations of independent Poisson variates $Y_{k} \sim Poi (μ_{k})$ with $μ_{k} = r_{k} e^{η_{k}}$ , where $η_{k} = \ln λ_{k}$ (Holford 1980; Laird and Olivier 1981). This correspondence links hazard estimation to generalized linear models and their extensions.

Bin-specific hazard levels $λ_{k}$ do not represent a smooth hazard, but alternatively we can express the log hazard $η (t) = \ln λ (t)$ as a linear combination of smooth basis functions $B_{l} (t)$ : $η (t) = \sum_{l = 1}^{L} B_{l} (t) α_{l}$ . Evaluating $η ()$ at the bin midpoints $m_{k}$ leads to $η_{k} = \sum_{l = 1}^{L} B_{l} (m_{k}) α_{l} = \sum_{l = 1}^{L} b_{k l} α_{l}$ with known $b_{k l} = B_{l} (m_{k})$ . In this way, hazard estimation becomes Poisson regression, in which we estimate the regression coefficients $α_{l}$ . Eilers (1998) suggested using $B$ -splines as basis functions $B_{l} (t)$ and penalize the coefficients by a difference penalty to achieve proper smoothness. Extensions were presented in van Houwelingen and Eilers (2000), and a Bayesian formulation was presented in Lambert and Eilers (2005). We provide a more detailed account of the approach in the Appendix (online supplement).

The concepts of univariate hazard smoothing with $P$ -splines can be readily extended to two-dimensional hazard smoothing. We will apply these extensions to the estimation of $\overset{˘}{λ} (u, s)$ . If $λ (t, s)$ is required, then back-transformation via (2) will provide the estimate.

The basic ingredients stay the same: the data, now extending over two time axes, are binned into bins of equal size (now squares or rectangles), and the log-hazard, which is a bivariate surface, is represented via linear combinations of $B$ -splines so that the link to Poisson regression persists. Difference penalties on the regression parameters again control the smoothness of the estimated hazard.

First, the two axes $u$ and $s$ are divided into $n_{u}$ and $n_{s}$ bins of equal length, which extend over the range of observed values so that, the $(u, s)$ -plane is covered by an $n_{u} \times n_{s}$ tessellation of small squares (or rectangles if bin width differs between $u$ - and $s$ -axis), which in the $(t, s)$ -plane correspond to parallelograms. The information of an individual can now be represented by two $n_{u} \times n_{s}$ matrices, the first one containing the time at risk in each square, and the second containing the information on the event. If the individual experiences an event, then the cell where the event is located is equal to 1. If the individual is right-censored, all cells are zero. From these individual matrices, the aggregated matrices of total exposure times $R = [r_{j k}]$ and total event counts $Y = [y_{j k}]$ are calculated.

The values $η_{j k} = \log {\overset{˘}{λ}}_{j k}$ are the log-hazard evaluated at the centers of the two-dimensional bins. These midpoints $m_{j k}$ have two components $m_{j k} = (m_{j}^{u}, m_{k}^{s})$ one coordinate for each axis. The matrix $E = [η_{j k}] = [\log {\overset{˘}{λ}}_{j k}]$ is of dimension $n_{u} \times n_{s}$ , too.

The log-likelihood is equivalent to that of a Poisson model for the counts $y_{j k}$ , where:

Y \sim Poisson (M) M = [μ_{jk}] = R ⊙ [{\overset{˘}{λ}}_{jk}] .

(3)

Here,

\overset{˘}{λ}

indicates that we are working in the transformed

(u, s)

-plane and

⊙

indicates element-wise multiplication.

The log-additive model (A) and the varying-coefficient model (B), see Section “Temporal Change: Time Scales, Covariates, Effects,” are straightforward extensions of the univariate $P$ -spline regression described in the Appendix (online supplement) , see Eilers and Marx (2021: Chapter 4), .

The log-additive model (A)

η (u, s) = θ (s) + ϕ (u)

(4)

requires two smooth but univariate components. The functions $ϕ (u)$ and $θ (s)$ are specified as linear combinations of two marginal $B$ -spline bases

θ (s) = \sum_{j = 1}^{c_{s}} B_{k}^{s} (s) α_{j}^{s} a n d ϕ (u) = \sum_{j = 1}^{c_{u}} B_{j}^{u} (u) α_{j}^{u} .

(5)

The data are vectorized, that is,

y = vec [Y], r = vec [R]

and

η = vec [E],

and we can express the linear predictor as

η = B α

, where the

n_{u} n_{s} \times (c_{u} + c_{s})

matrix

B

B = [B_{u} | B_{s}]

and the

(c_{u} + c_{s}) \times 1

vector

α

α = [\begin{matrix} α^{u} \\ α^{s} \end{matrix}] .

The rows of $B_{u}$ and of $B_{s}$ contain the marginal basis functions evaluated at the points $m_{j}^{u}$ and $m_{k}^{s}$ that correspond to the position of $y_{j k}$ in the respective row. The two components of $α$ are penalized separately with different smoothing parameters $ϱ_{u}$ and $ϱ_{s}$ . The $(c_{u} + c_{s}) \times (c_{u} + c_{s})$ penalty matrix $P$ is a blockdigonal matrix $P = diag (P_{u}, P_{s})$ with $P_{u} = ϱ_{u} D_{u}^{T} D_{u}$ and $P_{s} = ϱ_{s} D_{s}^{T} D_{s}$ . The components of an additive model are defined only up to a constant, which leads to a singular system (introducing an identifiability problem). Adding a small ridge penalty and using $P^{+} = P + κ I$ as penalty matrix removes this identifiability problem. The value $κ$ is small, typically in the order of $κ = 10^{- 6}$ .

With this particular structure, we can solve the system of equations obtained from maximizing the penalized log-likelihood of the model (as derived in Appendix equation (A.9), online supplement), this time depending on two smoothing parameters $ϱ_{u}$ and $ϱ_{s}$ . The partitioned structure of the model also allows partitioning of the effective dimension to give partial EDs for each component. For details, see Eilers and Marx (2021).

For the varying coefficient model (B)

η (u, s) = θ (s) + ϕ (s) u,

(6)

η (u, s) = θ (s) + ϕ (s) u

, only slight modifications are needed. A single marginal

B

-spline basis (along

s

) is required and

θ (s)

is identical to the representation in (5). To not confuse coefficient vectors, we denote the vector by

α^{θ}

(instead of

α^{s}

). The component

ϕ (s) u

can also be expressed using the same basis, however, with a different coefficient vector

α^{ϕ}

. The rows of the resulting matrix

B_{s}

for this component have to be multiplied by the

u

-values that correspond to the position of the respective element

y_{j k}

, that is,

m_{j}^{u}

. If

U

denotes an

n_{u} n_{s} \times n_{u} n_{s}

diagonal matrix that has these values in the corresponding rows, we obtain for the varying coefficient model a regression matrix

B = [B_{s} | U B_{s}]

and a coefficient vector

α = [\begin{matrix} α^{θ} \\ α^{ϕ} \end{matrix}]

of length

2 c_{s}

. The penalty matrix again is blockdiagonal:

P = diag (ϱ_{θ} D_{s}^{T} D_{s}, ϱ_{ϕ} D_{s}^{T} D_{s})

. And again, we arrive at a system of the form (A.9) (Appendix, online supplement), governed by two smoothing parameters

ϱ_{θ}

and

ϱ_{ϕ}

The most flexible model (C) is a smooth bivariate surface $η (u, s)$ . This is obtained by using tensor products of $B$ -splines as a basis. Tensor products of $B$ -splines extend the idea illustrated in Appendix Figure A1 (online supplement), to two dimensions. In Figure 2, a small subset of such a tensor product basis is shown. Each two-dimensional basis function will be multiplied by a coefficient; the number of basis functions can be relatively large, but the coefficients will be penalized by difference penalties to prevent overfitting.

Figure 2.

An example of a $3 \times 3$ tensor product basis of $B$ -splines. The tensor product basis contains a total of 9 $B$ -splines. In this illustration, the basis functions are well separated to show the principle. In practice, the basis functions are denser and overlap to ensure smooth results, see the one-dimensional case in Appendix Figure A1 (online supplement).

More specifically, the model is as follows. Two marginal $B$ -spline bases are constructed (see Appendix equation (A.5), online supplement), one for each of the two time dimensions $u$ and $s$ . They are denoted by $B_{u}$ and $B_{s}$ , and their dimension is $n_{u} \times c_{u}$ and $n_{s} \times c_{s}$ , respectively.

Their tensor product

B = B_{s} \otimes B_{u}

(7)

is of dimension $n_{u} n_{s} \times c_{u} c_{s}$ , the symbol $\otimes$ denotes the Kronecker product. The matrix $B$ , which is bi-dimensional and has full local support (Currie et al. 2004), is used as a regressor matrix in the Poisson model (3). For the additive model (4), $c_{u} + c_{s}$ coefficients are needed, the interaction model has $c_{u} c_{s}$ coefficients.

In this way, the model for the log-hazard becomes $η = B α$ , where $η = vec (E)$ . The vector operator stacks the columns of a matrix into a vector. Thus, the use of tensor product $B$ -splines and the vector operator leads us to the same equations as for the $P$ -spline model for one-dimensional hazards presented in the Appendix, online supplement.

While the extension to two dimensions is relatively straightforward, the notation can become cumbersome, and a clever arrangement of quantities will allow efficient computations. The $c_{u} c_{s}$ coefficients for the tensor matrix $B$ can be organized in a matrix $A = [α_{l m}]$ of dimensions $c_{u} \times c_{s}$ . The linear predictor $η_{j k} = \log λ_{j k}$ can then be expressed as a linear combination of the tensor product basis $B$ and the $c_{u} c_{s}$ vector of associated coefficients $α$ as:

η_{j k} = \log {\overset{˘}{λ}}_{j k} = \sum_{l = 1}^{c_{u}} \sum_{m = 1}^{c_{s}} b_{j l}^{u} b_{k m}^{s} α_{l m},

(8)

which in matrix notation simply is $E = B_{u} A B_{s}^{T}$ , and which conveniently has the same dimension as the data matrices $R$ and $Y$ , that is, $n_{u} \times n_{s}$ .

In two dimensions, a penalty on the coefficients in $A$ is introduced. The columns and the rows of $A$ are penalized separately, by using

\,pen = ϱ_{u} \sum_{l = 1}^{c_{u}} ‖ D_{u} α_{: l} ‖^{2} + ϱ_{s} \sum_{m = 1}^{c_{s}} ‖ D_{s} α_{m :}^{T} ‖^{2} = ϱ_{u} \sum_{l = 1}^{c_{u}} α_{: l}^{T} D_{u}^{T} D_{u} α_{: l} + ϱ_{s} \sum_{m = 1}^{c_{s}} α_{m :} D_{s}^{T} D_{s} α_{m :}^{T},

(9)

where the first term in the sum penalizes the columns of $A$ and the second term penalizes the rows of $A$ (Currie et al. 2004). Here, $α_{: l}$ indicates column $l$ of $A$ and $α_{m :}$ indicates row $m$ of $A$ . The matrices $D_{u}$ and $D_{s}$ are difference matrices of order $d$ (see also Appendix equation (A.7), online supplement), referring to the $u$ (age at entry) and $s$ (duration) axes. Finally, $ϱ_{u}$ and $ϱ_{s}$ are the smoothing parameters which can have different values to allow different amounts of smoothing along the rows and columns of $A$ (referred to as “anisotropic smoothing”). In matrix notation, equation (9) translates to $P = ϱ_{u} I_{s} \otimes D_{u}^{T} D_{u} + ϱ_{s} D_{s}^{T} D_{s} \otimes I_{u}$ with $I_{u}$ and $I_{s}$ being identity matrices of appropriate dimensions.

The optimal values of the smoothing parameters are found by minimizing the Akaike information criterion (AIC) (or the Bayesian information criterion (BIC)), as a function of the pair ( $ϱ_{u}$ , $ϱ_{s}$ ). This is achieved either by grid search or by numerical minimization. After the optimal values for $ϱ_{u}$ and $ϱ_{s}$ are selected, the estimated regression coefficients are inserted in equation (8), and $\overset{˘}{λ} (u, s)$ is estimated by:

\hat{\overset{˘}{λ}} (u, s) = \exp [B_{u} \hat{A} B_{s}^{T}] .

(10)

Uncertainty estimates are derived as byproducts of the iterative weighted least square algorithm. The variance–covariance matrix of $\hat{α}$ is

Cov (\hat{α}) \approx (B^{T} \hat{W} B + P)^{- 1},

(11)

where $B$ is the tensor product basis, $\hat{W}$ is the matrix containing the weights, and $P$ the penalty matrix. The variance–covariance matrix of the estimated log-hazard is $Cov (\hat{η}) = Cov (B \hat{α}) = B Cov (\hat{α}) B^{T}$ , the square root of its diagonal elements give the standard errors of the $\hat{η}$ . Standard errors of the hazard $\hat{λ}$ are obtained by the delta method: $s . e . (\hat{λ}) = e^{\hat{η}} s . e . (\hat{η})$ . The asymmetric confidence interval for the hazard $\hat{λ}$ can be obtained as $\exp (\hat{η} \pm 1.96 \times s . e . (\hat{η}))$ .

Transitions From Cohabitation to Marriage or Separation

Background

Premarital cohabitation, that is, living together with a romantic partner without being married, is a common form of living arrangement. In many countries in Europe and in the United States, premarital cohabitation has replaced direct marriage in the process of union formation (Sassler and Lichter 2020; Sobotka and Berghammer 2021). Therefore, scholars are increasingly interested in understanding the timing and determinants of transitions out of unmarried cohabiting unions. The majority of cohabiting individuals will eventually marry their partner (Nazio and Blossfeld 2003; Sassler and Lichter 2020), while some will experience the dissolution of their cohabitation (Hiekel et al. 2014). In order to provide a comprehensive picture of these processes, researchers need to properly account for two time scales, which are drivers of such transitions: duration of the cohabitation and age of the individual. However, published studies account for the presence of these two time scales only partly.

Research focused on the transitions out of cohabitation, either into marriage or into union dissolution, usually considers union duration as the most important time scale (Lu et al. 2012), while controlling for age at union formation (Hiekel et al. 2015; Wright 2019) or both age and age at start of cohabitation (Le 2002). Union duration is often considered as the main time scale for several reasons. Unmarried cohabitation can be seen as a trial marriage, in which individuals test the suitability of the partner as a spouse (Klijzing 1992). Since cohabiting couples learn about their compatibility over time, longer cohabitations might increase the risk of transitioning to marriage (Schnor 2014), while the risk of separation is higher at shorter durations. The latter is especially true for people who have started a cohabitation at a younger age, and therefore might be more likely to be in a poor match (Ermisch and Francesconi 2000).

Chronological age might be a proxy for internalized or ascribed norms for behaviors regarding life transitions (Settersten and Mayer 1997). Individuals often organize their lives around chronological age-deadlines (Settersten and Hagestad 1996). As an example, in a social context where childbearing is more widely accepted inside a marriage, such as in Western Germany, individuals might feel pressured to marry before they reach a certain age, to then transition to parenthood. Additionally, younger couples might be less likely to transition to marriage and more likely to undergo a union dissolution (Ermisch and Francesconi 2000). Hence, many authors study life transitions using chronological age as the time scale.

Following these theoretical considerations, we argue that both age and union duration are important time scales that influence the individuals’ decision to marry or to dissolve their cohabitation. Moreover, the hazard of these life events might be quite different for different combinations of age and duration of the cohabitation, indicating possibly complex interactions between the two time dimensions. For example, individuals who enter a cohabitation with their partner at a younger age might be less inclined to marry within a couple of years and therefore remain in unmarried cohabitation for longer times. Contrary to this, individuals and especially women, who start cohabiting in their thirties, might feel more (biological or societal) pressure to marry their partner quickly to successively transition to childbearing.

In the following, we show how models (A) to (C) described in Section “Temporal Change: Time Scales, Covariates, Effects” capture the interplay of age and duration in the hazard of ending an unmarried cohabitation, in either marriage or separation, and which level of complexity provides the best fit.

Data and Study Design

We use data from waves 1 to 11 of the pairfam (Panel Analysis of Intimate Relationships and Family Dynamics), release 11.0 (Brüderl et al. 2020a). A detailed description of the study can be found in Huinink et al. (2011).

Pairfam is a longitudinal panel survey, conducted annually, providing rich data on the formation and development of intimate relationships and families. The panel started with about 12,000 randomly selected individuals (anchors) of three different birth cohorts (1991–1993, 1981–1983, and 1971–1973) living in West and East Germany. The first wave of interviews was conducted in 2008. In 2009, about 1,500 individuals living in East Germany were sampled as part of the panel study DemoDiff, which was initiated following the design of pairfam. Beginning with wave five, the two studies have been fully integrated under pairfam. Finally, in 2019 (wave 11), a refreshment sample was drawn, which added new individuals from cohorts 1981–1983 and 1991–1993 together with a new cohort of individuals born in 2001–2003.

In this study, we use the dataset biopart, which is a ready-to-use dataset prepared by the pairfam team. The dataset biopart provides both retrospective and prospective information on individuals’ relationship histories from age 14, including cohabitations and marriages, on a monthly basis, and it is constructed with data collected in each wave. Details of biopart, as well as of the study design, can be found in the Data Manual (Brüderl et al. 2020b).

For this analysis, we included men and women living in East or West Germany who had experienced at least one non-marital cohabitation with a romantic partner of the opposite sex during the study period. At their first interview, individuals are asked to list all romantic relationships from age 14, together with the dates of important events (beginning of cohabitation, marriage, union dissolution, etc.). Therefore, information on all cohabiting unions before and after the first interview is available in the biopart dataset.

In our analysis, we include two sets of cohabitations: First, we include all unmarried cohabiting unions that were formed before but have not ended by the time of the first interview. For these cohabiting unions, the date of the first interview is taken as entry into the study. The cohabitation has been running at that point for a duration $\underline{s} > 0$ , and the observation is left truncated. This implies that we exclude all cohabiting unions that have started and ended before the individual was first interviewed. Second, we also include all cohabitations that have started some time after the individual was interviewed for the first time. In these cases, the observation is not left truncated, as we are able to follow the cohabitation from the start, and hence $\underline{s} = 0$ . Finally, we consider only one cohabitation per individual in the sample.

These choices are visualized in Figure 3: here we represent all transitions between states that might be observed thanks to the design of pairfam and the data available in the biopart dataset.

Figure 3.

Sketch of biopart trajectories that were included in the analyses. C = Start of cohabitation, M = Marriage, S = Separation. Cohabiting unions that ended before the first interview were excluded. Cohabitations that started before the first interview were included and delayed entry (left truncation) was accounted for. Cohabitations that started during the observation period were also included in the analysis.

In the biopart dataset, there are 14,235 cohabitations, of which 8,451 had already ended before the first interview. Same-sex couples were excluded because regulations regarding same-sex marriages have changed during the study period (172 cohabitation). Additionally, we excluded 58 cohabitations that ended because of the partner’s death and 62 cohabitations with unclear ending dates. We also excluded cohabitations that started directly with a marriage (143 cohabitations). Finally, we only considered cohabitations that started when the individual was 18 years or older, which defines full legal age in Germany. Because of this age restriction, only one individual from the youngest cohort (2001–2003) would have entered our sample, so we decided to restrict to the three oldest cohorts.

The key variables in our analysis are date of birth, date of entry into cohabitation, date of marriage, date of union dissolution, dates of first and last interviews, sex, and region of residence at each interview. The date variables are used to identify ages at each transition and durations of the cohabitation. Region of residence is collected at each wave, and distinguishes between East and West Germany. This variable takes into account that individuals can move between the two regions, and can be different from the region of birth. However, we consider the region of residence at entry into cohabitation, and keep it as a fixed covariate.

We restrict the sample to individuals with non-missing values for all of the covariates (removing 16 records). All individuals are censored at the time of their last interview, if they did not experience any event before then. The maximum observed duration of the cohabitation was a bit less than 27 years; however, events are extremely rare after 20 years of cohabitation, and only very few individuals are still at risk (Supplementary Figure S2, online supplement). We censored three individuals, two women from Eastern Germany and one woman from Western Germany, who had events after more than 23 years of cohabitation, as they were considered outliers with respect to the key variable in our analysis (duration of the cohabitation), therefore removing two marriages and one separation.

Finally, we only consider one cohabitation per individual, which is the first cohabitation in the data that matches the criteria described above, arriving at a final sample size of 4,669 observations. The data selection process is also described in a flowchart that we include in Supplementary Figure S1 (online supplement).

In all our analyses, we consider four different groups defined by the sex and the region of residence of the respondent. Table 2 provides some basic information.

Table 2.

Descriptive Statistics of the Analyzed Sample.

Sample	Women West	Women East	Men West	Men East
$N$ (at baseline)	1,649	852	1,443	725
Age at cohabitation	26.2 (5.8)	26.4 (6.1)	27.3 (5.6)	27.9 (5.6)
Married (%)	372 (22.6)	216 (25.4)	364 (25.2)	172 (23.7)
Age at marriage	31.1 (5.5)	32.4 (5.6)	32.4 (5.6)	33.8 (5.3)
Duration at marriage	4.3 (3.5)	5.7 (3.9)	4.2 (4.0)	5.7 (4.2)
Separated (%)	324 (19.7)	169 (19.8)	276 (19.1)	131 (18.1)
Age at separation	27.9 (6.8)	29.4 (7.4)	28.3 (6.2)	31.7 (6.3)
Duration at separation	2.8 (2.8)	4.0 (4.2)	2.3 (2.3)	3.6 (3.5)
Censored (%)	953 (57.8)	467 (54.8)	803 (55.7)	422 (58.2)

Note. First row presents the sample sizes of the four groups (at entry into the risk set). The other numbers are mean and standard deviation (in parentheses). Age at marriage/separation and duration of cohabitation at marriage/separation are calculated by restricting the sample to individuals who have experienced the respective event (excluding ongoing cohabitations).

Estimated Hazards Over two Time Scales

We present the results separately for the event “marriage” and for the event “separation,” and we compare the interaction model (8) with the simpler additive model (4) and the varying-coefficient model (6). In the analysis, $t$ is the age of the individual, $s$ is the duration of the cohabitation, and $u = t - s$ is the age at the beginning of the cohabitation. We measure all variables in years. All models are estimated over the $(u, s)$ -plane (as described in Section “Smoothing Bivariate Hazard Rates with P-splines”). If the competing event is experienced, then an observation is censored. Censoring is also introduced if the individual was still cohabiting at the last available interview date. Censoring times are, in this case, calculated from the date of the last interview.

We estimate the smooth hazard of marriage or separation for each of the four groups (men/women; East/West) separately. In all models, we choose cubic $B$ -splines (degree $p = 3$ ) and second-order difference penalties ( $d = 2$ ). The data are binned in intervals of length 0.5 (years). Optimal values of the smoothing parameters $ϱ_{u}$ and $ϱ_{s}$ are obtained by numerical minimization of the AIC using function ucminf in R-package ucminf (Nielsen and Mortensen 2022).

Transitions From Cohabitation to Marriage

Table 3 summarizes the optimal values of the smoothing parameters, which were obtained by minimizing AIC, the resulting effective dimensions, and the values of AIC for the twelve models. In each case, the smallest AIC is marked in bold.

Table 3.

Hazard of Marriage.

	Log-additive					Varying-coefficient					Full bivariate
	$\lg ϱ_{θ}$	$\lg ϱ_{ϕ}$	${ED}_{θ}$	${ED}_{ϕ}$	AIC	$\lg ϱ_{θ}$	$\lg ϱ_{ϕ}$	${ED}_{θ}$	${ED}_{ϕ}$	AIC	$\lg ϱ_{u}$	$\lg ϱ_{s}$	$ED$	AIC
MW	1.0	1.0	4.7	3.5	853	5.0	3.8	2.0	4.0	840	0.2	$- 0.3$	15.2	808
ME	3.0	1.5	1.7	2.8	672	1.8	5.0	2.8	2.3	669	1.7	0.4	6.9	670
WW	$- 0.5$	1.0	8.1	3.7	865	0.8	5.5	5.2	2.0	858	0.4	0.1	12.9	827
WE	3.0	$- 0.5$	1.7	6.1	772	3.5	3.3	2.0	5.8	766	1.7	$- 0.1$	8.0	767

Note. Optimal smoothing parameters ( $\lg = \log_{10}$ ), resulting EDs, and AIC for the log-additive, varying-coefficient, and full bivariate model. MW = Men West, ME = Men East, WW = Women West, WE = Women East; ED=effective dimensions; AIC = Akaike information criterion.

In all four groups, a log-additive (proportional hazards) model is too simple. A single hazard shape over $s$ is not sufficient to capture the interplay between age and duration of cohabitation. When comparing the varying-coefficient model and the full bivariate model, we see that both for men and for women in West Germany, the latter clearly wins out. In East Germany, the varying-coefficient model and the full bivariate model fare almost equally. So in East Germany, we would select the varying-coefficient model as final model.

Figure 4 shows the estimated surfaces of $\overset{˘}{λ} (u, s)$ for men and the three different models. In Figure 5, the hazard $\overset{˘}{λ} (u, s)$ of marriage over duration $s$ is shown for selected values of $u$ , the age at start of cohabitation ( $u = 20, 25, \dots, 45$ ). The layout matches the one in Figure 4. The corresponding plots for women, which lead to rather similar conclusions, are deferred to Supplementary Figure S5 (online supplement).

Figure 4.

Hazard of marriage over age at entry into and duration of cohabitation: Comparison of the log-additive model (left), varying-coefficient model (center), and 2D interaction model (right) for Men West (top) and Men East (bottom). Corresponding values of effective dimensions and AIC are in Table 3.

Figure 5.

Hazard of marriage over duration of cohabitation for selected values of age at entry into cohabitation ( $u$ ) for men. The layout corresponds to Figure 4. (Colour coding of the hazard lines is given in the bottom left panel and applies to all subplots.)

In general, hazard levels for marriage are lower in East Germany than in West Germany, confirming previous research. For West German men, there is a noticeable increase in the hazard for cohabitations that start between ages 25 and 30 and end, due to marriage, after about five years (see top right panel). For cohabitations that are entered at later ages, marriage rates also peak at about four to five years, although at lower levels. Also, there is an increase in the marriage hazard for individuals who started cohabiting quite young (and did not split up) with increasing duration of their cohabitation, peaking at about 10–15 years, for $u = 20$ , pointing towards a link with family formation. The varying-coefficient model (top row, center) attempts to capture these features; however, its interaction structure is not sufficiently flexible to catch the age-duration dynamics that require a shift in the mode as well as changing levels and varying declines of marriage hazards after the mode (for men in West Germany). As a consequence, when applied to these data, the varying-coefficient model overestimates the hazard of marriage for ages $u$ above 35 years.

For men in East Germany, the plots in the bottom row, center, and right, of Figures 4 and 5 confirm the results in Table 3: The estimates of $\overset{˘}{λ} (u, s)$ are quite similar, so in this case the simpler interaction structure of the varying-coefficient model is sufficiently flexible.

The estimated surfaces $λ (t, s)$ over currently attained age $t$ and duration of cohabitation $s$ , which are obtained by using equivalence (2), are given in Supplementary Figures S3 and S4 (online supplement).

Transitions From Cohabitation to Separation

The results for the hazard of separation in the four groups are summarized in Table 4.

Table 4.

Hazard of Separation.

	log-additive					varying-coefficient					full bivariate
	$\lg ϱ_{θ}$	$\lg ϱ_{ϕ}$	${ED}_{θ}$	${ED}_{ϕ}$	AIC	$\lg ϱ_{θ}$	$\lg ϱ_{ϕ}$	${ED}_{θ}$	${ED}_{ϕ}$	AIC	$\lg ϱ_{u}$	$\lg ϱ_{s}$	$ED$	AIC
MW	2.3	5.0	2.4	1.5	556	4.0	6.5	2.0	2.0	559	1.5	5.1	5.0	557
ME	3.2	5.0	1.6	1.5	528	3.8	5.5	2.0	2.2	531	0.9	3.0	5.4	530
WW	0.6	2.6	5.4	1.7	668	5.5	5.0	2.0	2.6	690	0.1	0.1	12.8	669
WE	1.5	2.4	3.2	1.8	576	3.3	5.3	2.1	2.5	589	0.4	2.1	6.8	578

Different from the marriage rates, the conclusions here are uniform and lead to a much simpler result: In all four groups, the additive model has the smallest AIC, so a single hazard $λ_{0} (s) = \exp {θ (s)}$ changes proportionally with $u$ .

Figure 6 presents the results for men. The estimated functions $θ (s)$ are linear, confirmed by the central column in the figure, so the hazard to separate declines exponentially with duration $s$ of the cohabitation—a simple parametric model will do. The more flexible bivariate hazard model (right column) basically results in the same hazard shapes, without strict proportionality, though. The hazard of separation in general is highest when cohabitation is entered young and declines with increasing age at the start of cohabitation. Additional illustrations, such as the surfaces $λ (t, s)$ , and the corresponding results for women can be found in the online supplement.

Figure 6.

Hazard of separation over duration of cohabitation for selected values of age at entry into cohabitation ( $u$ ), men in West Germany (top row) and East Germany (bottom row). Left: Hazard of the log-additive model (minimal AIC, cf. Table 4). Center: Log-hazard confirming the exponential shape of the hazard. Right: Hazard of the full bivariate model. Note. AIC=Akaike information criterion.

In summary, the proposed methodology identifies quite different models for the hazard of marriage and of separation. The level of complexity ranges from a simple PH model to a full bivariate hazard surface. AIC (or BIC) provides a convenient measure for model comparison and selection, effective dimensions of the models give further useful information on the model complexity. When a simpler model is selected as optimal, then the full bivariate surface reproduces the simpler hazard.

Discussion and Outlook

When processes evolve in time, several time scales may be important, and models that are used to analyze such processes should allow us to incorporate more than one time scale. This will allow us to investigate their joint roles in the process. This paper offers an approach to incorporate multiple time scales in the analysis of time-to-event data and demonstrates its use in the study of transitions out of nonmarital cohabitation to either marriage or separation. We review how the simultaneous analysis of two time scales in hazard models can be reformulated using an “age-period-cohort”-like relationship. The impact of two time scales can then be expressed in a single time scale model with covariate effects of different complexity, ranging from a proportional hazards structure to a full bivariate hazard surface in which the covariate modulates the hazard smoothly but without additional imposed structure.

To estimate these hazard models, we suggest to use $P$ -splines. The data are binned on a fine grid and the different hazard models are represented, on log-scale, by linear combinations of one- or two-dimensional $B$ -spline bases. Difference penalties on the coefficient vectors control the smoothness of the hazard surface. This approach links hazard estimation to Poisson regression for which efficient fitting algorithms exist. Effective dimensions of the fitted models are automatically provided and allow us to compare models of different complexity via AIC (or BIC) in a straightforward manner.

In the application, we demonstrate how the proposed approach identifies, using a single methodology, interactions of different complexity (for the hazard of marriage) or much simpler proportional hazards (for transitions to separation), which can even be modeled parametrically.

The initial binning of the data is sometimes claimed to be a drawback even though it is this binning that was recognized, already several decades ago, to provide the extremely productive link between hazard modeling and Poisson regression. Even somewhat narrow bins that contribute to the flexibility of $P$ -spline models do not pose obstacles regarding data size or computation times, given current technology.

Recent alternative approaches to model hazards over two time scales use flexible parametric models (Batyrbekova et al. 2022; Bower et al. 2022). Full bivariate hazard surfaces, which are required for some hazards in our application, are mentioned but, as far as we can tell, have not been implemented yet.

Our analysis of transitions from cohabitation to marriage or separation in Germany shows interesting results. Previous research, using the same data, found that the rate of marrying (or separating) steadily increases in the first five years after cohabitation and thereafter remains unchanged, and that age at cohabitation did not significantly change the hazard rates (Hiekel et al. 2015). Our results tell a more complicated story.

In this study, the focus was on assessing the impact of several time scales, and no further covariates were included. It is, however, possible to use hazards over two time scales as a baseline hazard in a proportional hazards (PH) model. The extension to the PH regression model is discussed in Carollo et al. (2025), and it is implemented in the R-package, currently only for fixed covariates. We decided to omit details about the implementation of the PH model; however, we do present the results of PH analysis in the online supplement, using the same data as in Hiekel et al. (2015), including the same covariates but using a two-dimensional baseline hazard instead. We find that, even when using the same sample and controlling for the same covariates, the baseline hazard of marriage and of separation over the duration of the cohabitation differ substantially by age at entry into cohabitation. This comparison provides evidence that differences between the two analyses may be explained by the use of the considerably more flexible hazard models, rather than differences in the analyzed sample.

Our results also show differences in the rate of marriage between East and West German men and women, which are coherent with the literature on the topic. Klarner and Knabe (2017) found that East German individuals see marriage as an institution with little relevance for everyday behavior or for the meaning of their partnerships, while Western German couples still believe in marriage as an institution, and one of the biggest arguments in favor of marriage is the formalization of the relationship.

We analyzed cohabitation from age 18 onward. To check how sensitive our results are to this choice, we repeated the analysis, enlarging our sample to ages at entry into cohabitation as low as 15 years by adding 230 cohabitations. The patterns shown in the hazard of marriage or separation remain mostly unchanged. Also, we considered one cohabitation per individual in our sample, which might not have been the first ever cohabitation experienced by the individual. It is possible that the risk of marriage or separation among first-time cohabiters is different from those who are in their second or higher cohabitation. To check if our results are robust regarding the order of the cohabitation, we ran another analysis restricting the sample to the first-ever cohabitation only, thereby removing about 25.5% of all observations. We again find that the way both time dimensions interact is the same among first-time cohabiters and among all cohabiters. Both these additional sensitivity analyses are presented in the online supplement.

The analyses presented in this paper were produced using the TwoTimeScales R-package (Carollo et al. 2024), available on CRAN. The pairfam data are, in principle, available to the interested reader, but have to be applied for (see https://www.pairfam.de/en). Similar analyses, using publicly available data, are provided in vignettes to the R package.

The current model is limited to the inclusion of two time scales. Occasionally, additional time scales might be relevant, although we suspect that more than three time scales may challenge our power of imagination. Extension of the $P$ -spline models to more than two dimensions is in principle possible and has been used in the literature (see Currie et al. 2006). Extending the present approach to more than two time scales will, however, require additional work regarding computations and presentation of results. Likewise, extending the methodology to regression models that include additional covariates other than in a linear predictor is considered future research.

We analyzed cause-specific hazards, considering the competing events of marriage and separation. From the cause-specific hazards, transition probabilities can be calculated, as shown in Carollo et al. (2024), to complete a competing risks analysis. Extensions to more complex multi-state models will be the topic of future research, too. In view of the current debate about the role that childbearing has in union formation (see, e.g., Wright 2019, Schnor 2014, and Lappegård et al. 2018) such extensions are of practical relevance. It is conceivable that with additional transitions, also additional time scales gain importance.

Supplemental Material

sj-pdf-1-smr-10.1177_00491241251374193 - Supplemental material for Analysis of Time-to-Event Data With Two Time Scales. An Application to Transitions out of Cohabitation

Supplemental material, sj-pdf-1-smr-10.1177_00491241251374193 for Analysis of Time-to-Event Data With Two Time Scales. An Application to Transitions out of Cohabitation by Angela Carollo, Hein Putter, Paul H.C. Eilers and Jutta Gampe in Sociological Methods & Research

Supplemental Material

sj-pdf-2-smr-10.1177_00491241251374193 - Supplemental material for Analysis of Time-to-Event Data With Two Time Scales. An Application to Transitions out of Cohabitation

Supplemental material, sj-pdf-2-smr-10.1177_00491241251374193 for Analysis of Time-to-Event Data With Two Time Scales. An Application to Transitions out of Cohabitation by Angela Carollo, Hein Putter, Paul H.C. Eilers and Jutta Gampe in Sociological Methods & Research

Footnotes

Acknowledgments

We are very grateful to Nicole Hiekel from the Max Planck Institute for Demographic Research for her helpful comments on an earlier version of this manuscript. Her insights greatly helped to shape the presentation of the application and discussion of the findings. This paper uses data from the German Family Panel pairfam, coordinated by Josef Brüderl, Sonja Drobnič, Karsten Hank, Bernhard Nauck, Franz Neyer, and Sabine Walper. pairfam is funded as a long-term project by the German Research Foundation (DFG).

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Preregistration Statement

This study was not preregistered.

Supplemental Material

Supplemental material and Appendix for this article are available online.

ORCID iDs

Angela Carollo

Hein Putter

Paul H.C. Eilers

Jutta Gampe

Data and Code Availability Statement

The data were used under license from GESIS for the current study and are not publicly available. However, the data can be accessed using the instructions provided here: https://www.pairfam.de/en/data/data-access/. The code used for this study and documentation for the code are available in the GitHub repository (Carollo, 2025) indexed with Zenodo¹.

Notes

Author Biographies

Angela Carollo is a Statistical Analyst in the laboratory of Fertility and Well-Being of the Max Planck Institute for Demographic Research (MPIDR) in Rostock (Germany).

Jutta Gampe is the Head of the Laboratory of Statistical Demography at the Max Planck Institute for Demographic Research (MPIDR) in Rostock (Germany).

Hein Putter is full professor at the Leiden University Medical Center in Leiden (The Netherlands). He is coauthor of the book “Dynamic Prediction in Clinical Survival Analysis” (CRC Press, 2011), with Hans van Houwelingen.

Paul H.C. Eilers is an emeritus professor at Erasmus University Medical Center in Rotterdam (The Netherlands). Together with Brian Marx, he wrote the book “Practical Smoothing. The Joys of P-splines” (Cambridge, 2021).

References

Batyrbekova

Bower

Dickman

P. W.

Ravn Landtblom

Hultcrantz

Szulkin

Lambert

P. C.

Andersson

T. M.-L.

. 2022. “Modelling Multiple Time-scales with Flexible Parametric Survival Models.” BMC Medical Research Methodology 22(1): 290–302.

Bender

Groll

Scheipl

. 2018. “A Generalized Additive Model Approach to Time-to-event Analysis.” Statistical Modelling 18(3-4): 299–321.

Bender

Scheipl

. 2018. “Pammtools: Piece-wise Exponential Additive Mixed Modeling tools.” arXiv:1806.01042 [stat] (https://arxiv.org/abs/1806.01042).

Blossfeld

H.-P.

Hamerle

Mayer

K. U.

. 1989. Event History Analysis. New York: Lawrence Erlbaum Associates.

Bower

Andersson

T. M-L.

Crowther

M. J.

Lambert

P. C.

. 2022. “Flexible Parametric Survival Analysis with Multiple Timescales. Estimation and Implementation Using stmt.” The Stata Journal 22(3): 679–701.

Brüderl

Drobnič

Hank

Neyer

F. J.

Walper

Alt

Bozoyan

, et al. 2020a. “The German Family Panel (pairfam).” GESIS Data Archive, Cologne. ZA5678 Data file Version 11.0.0.

Brüderl

Frister

Hajek

Herzig

Lenke

Schütze

Schumann

. 2020b. pairfam Data Manual, Release 11.0. Munich: LMU Munich.

Carollo

. 2025. AngelaCar/SMR-Cohabitation2TS: v1.0.0-SMR-Cohabitation2TS (v1.0.0). Zenodo. https://doi.org/10.5281/zenodo.16738184.

Carollo

Eilers

P. H. C.

Gampe

. 2024. TwoTimeScales: Analysis of Event Data with Two Time Scales. R package version 1.0.0 (https://CRAN.R-project.org/package=TwoTimeScales).

10.

Carollo

Eilers

Putter

Gampe

. 2025. “Smooth Hazards With Multiple Time Scales.” Statistics in Medicine 44: e10297.

11.

Carollo

Putter

Eilers

Gampe

. 2024. “Competing Risks Models with Two Time Scales.” arXiv:2408.15058 [stat.ME] (https://arxiv.org/abs/2408.15058).

12.

Currie

I. D.

Durban

Eilers

P. H.

. 2004. “Smoothing and Forecasting Mortality Rates.” Statistical Modelling 4(4): 279–298.

13.

Currie

I. D.

Durbán

Eilers

P. H. C.

. 2006. “Generalized Linear Array Models with Applications to Multidimensional Smoothing.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(2): 259–280.

14.

Di Giulio

Impicciatore

Sironi

. 2019. “The Changing Pattern of Cohabitation: A Sequence Analysis Approach.” Demographic Research 40: 1211–1248.

15.

Efron

2002. “The Two-way Proportional Hazards Model.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64(4): 899–909.

16.

Eilers

P. H. C.

1998. “Hazard Smoothing with

B

-splines.” Pp. 200-207 in Proceedings of the 13th International Workshop on Statistical Modelling, edited by B.D. Marx and H. Friedl. Baton Rouge: Louisiana State University.

17.

Eilers

P. H. C.

Marx

B. D.

. 1996. “Flexible Smoothing with B-Splines and Penalties.” Statistical Science 11(2): 89–102.

18.

Eilers

P. H. C.

Marx

B. D.

. 2021. Practical Smoothing: The Joys of P-Splines. Cambridge: Cambridge University Press.

19.

Ermisch

Francesconi

. 2000. “Cohabitation in Great Britain: Not for Long, But Here to Stay.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 163(2): 153–171.

20.

Hiekel

Liefbroer

A. C.

Poortman

A.-R.

. 2014. “Understanding Diversity in the Meaning of Cohabitation Across Europe.” European Journal of Population 30(4): 391–410.

21.

Hiekel

Liefbroer

A. C.

Poortman

A.-R.

. 2015. “Marriage and Separation Risks Among German Cohabiters: Differences Between Types of Cohabiter.” Population Studies 69(2): 237–251.

22.

Hobcraft

Murphy

. 1986. “Demographic Event History Analysis: A Selective Review.” Population Index 52(1):3–27.

23.

Holford

T. R.

1980. “The Analysis of Rates and of Survivorship Using Log-linear Models.” Biometrics 36(2):299–305.

24.

Huinink

Brüderl

Nauck

Walper

Castiglioni

Feldhaus

. 2011. “Panel Analysis of Intimate Relationships and Family Dynamics (pairfam): Conceptual Framework and Design.” Zeitschrift für Familienforschung – Journal of Family Research 23: 77–101.

25.

Iacobelli

Carstensen

. 2013. “Multiple Time Scales in Multi-state Models.” Statistics in Medicine 32(30): 5315–5327.

26.

Keiding

1990. “Statistical Inference in the Lexis Diagram.” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 332(1627): 487–509.

27.

Klarner

Knabe

. 2017. “On the Normative Foundations of Marriage and Cohabitation: Results From Group Discussions in Eastern and Western Germany.” Demographic Research 36: 1637–1666.

28.

Klein

J. P.

Moeschberger

M. L.

. 2005. Survival Analysis. New York: Springer.

29.

Klijzing

1992. “‘Weeding’ in the Netherlands: First-Union Disruption Among Men and Women Born Between 1928 and 1965.” European Sociological Review 8(1): 53–70.

30.

Laird

Olivier

. 1981. “Covariance Analysis of Censored Survival Data Using Log-Linear Analysis Techniques.” Journal of the American Statistical Association 76(374): 231–240.

31.

Lambert

Eilers

P. H. C.

. 2005. “Bayesian Proportional Hazards Model with Time-varying Regression Coefficients: A Penalized Poisson Regression Approach.” Statistics in Medicine 24: 3977–3989.

32.

Lappegård

Klüsener

Vignoli

. 2018. “Why are Marriage and Family Formation Increasingly Disconnected Across Europe? A Multilevel Perspective on Existing Theories.” Population, Space and Place 24(2): e2088.

33.

Le Goff

J.-M

. 2002. “Cohabiting Unions in France and West Germany: Transitions to First Birth and First Marriage.” Demographic Research 7: 593–624.

34.

Lexis

1875. Einleitung in Die Theorie Der Bevölkerungsstatistik. Straßburg: De Gruyter.

35.

Qian

Cunningham

C.-L.

. 2012. “Estimating the Effect of Premarital Cohabitation on Timing of Marital Disruption: Using Propensity Score Matching in Event History Analysis.” Sociological Methods & Research 41(3): 440–466.

36.

Nazio

Blossfeld

H.-P.

. 2003. “The Diffusion of Cohabitation Among Young Women in West Germany, East Germany and Italy.” European Journal of Population – Revue Europeenne de Demographie 19(1): 47–82.

37.

Nielsen

H. B.

Mortensen

S. B.

. 2022. ucminf: General-Purpose Unconstrained Non-Linear Optimization. Version: 1.1-4.1.

38.

Page

H. J.

1977. “Patterns Underlying Fertility Schedules: A Decomposition by Both Age and Marriage Duration.” Population Studies 31(1): 85–106.

39.

Pencina

M. J.

Larson

M. G.

D’Agostino

R. B.

. 2007. “Choice of Time Scale and Its Effect on Significance of Predictors in Longitudinal Studies.” Statistics in Medicine 26(6): 1343–1359.

40.

Rebora

Galimberti

Valsecchi

M. G.

. 2015. “Using Multiple Timescale Models for the Evaluation of a Time-dependent Treatment.” Statistics in Medicine 34(28): 3648–3660.

41.

Sassler

Lichter

D. T.

. 2020. “Cohabitation and Marriage: Complexity and Diversity in Union-formation Patterns.” Journal of Marriage and Family 82(1): 35–61.

42.

Schnor

2014. “The Effect of Union Status At First Childbirth on Union Stability: Evidence From Eastern and Western Germany.” European Journal of Population 30(2): 129–160.

43.

Settersten

R. A.

. 2003. “Age Structuring and the Rhythm of the Life Course.” Pp. 81-98 in Handbook of the life Course. Springer.

44.

Settersten

R. A.

Hagestad

G. O.

. 1996. “What’s the Latest? Cultural Age Deadlines for Family Transitions.” The Gerontologist 36(2): 178–188.

45.

Settersten

R. A.

Mayer

K. U.

. 1997. “The Measurement of Age, Age Structuring, and the Life Course.” Annual Review of Sociology 23(1): 233–261.

46.

Sobotka

Berghammer

. 2021. Demography of Family Change in Europe. Cheltenham: Edward Elgar Publishing.

47.

Tassistro

Bernasconi

D. P.

Rebora

Valsecchi

M. G.

Anolini

. 2019. “Modeling the Hazard of Transition Into the Absorbing State in the Illness-death Model.” Biometrical Journal 62(3): 836–851.

48.

Thièbaut

A. C. M.

Bènichou

. 2004. “ Choice of Time-scale in Cox’s Model Analysis of Epidemiologic Cohort Data: a Simulation Study .” Statistics in Medicine 23(24): 3803–3820.

49.

van Houwelingen

H. C.

Eilers

P. H. C.

. 2000. “Non-proportional Hazards Models in Survival Analysis.” Pp. 151-160 in COMPSTAT, Physica-Verlag HD.

50.

Wolkewitz

Cooper

B. S.

Palomar-Martinez

Alvarez-Lerma

Olaechea-Astigarraga

Barnett

A. G.

Schumacher

. 2016. “Multiple Time Scales in Modeling the Incidence of Infections Acquired in Intensive Care Units.” BMC Medical Research Methodology 16(1): 116–126.

51.

Wright

2019. “Union Transitions and Fertility Within First Premarital Cohabitations in Canada: Diverging Patterns by Education?” Demography 56(1): 151–167.

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.35 MB

6.53 MB