Nonparametric double additive cure survival models: An application to the estimation of the non-linear effect of age at first parenthood on fertility progression

Abstract

This article introduces double additive models to describe the effect of continuous covariates in cure survival models, thereby relaxing the traditional linearity assumption in the two regression parts. This class of models extends the classical event history models when an unknown proportion of the population under study will never experience the event of interest. They are used on data from the German Socio-Economic Panel (GSOEP) to examine how age at first birth relates to the timing and quantum of fertility for given education levels of the respondents. It is shown that the conditional probability of having further children decreases with the mother's age at first birth. While the effect of age at first birth in the third birth's probability model is fairly linear, this is not the case for the second child with an accelerating decline detected for women that had their first kid beyond age 30.

Keywords

Bayesian P-splines births cure survival models continuous covariates double additive models fertility studies

1 Introduction

The increase in the age at first parenthood is a key indicator of the second demographic transition (Lesthaeghe, 1994, 2010; Sobotka, 2008; Sobotka et al., 2011). Albeit that the onset and pace of this development varied greatly across countries and in spite of the fact that some European countries have reported that postponement is even about to come to a halt, late first-time childbearing is a pertinent and common feature of the fertility patterns in European countries (Frejka and Sardon, 2006; Goldstein et al., 2009).

On the macro level, the great challenge has been to adequately account for the tempo effects that have distorted period fertility rates (Billari and Kohler, 2004; Sobotka, 2004). On the micro level, researchers have examined whether and how a late age at first birth affects parity progression, birth spacing and completed fertility(Marini and Hodsdon, 1981; Dommaraju, 2009; Bratti and Tatsiramos, 2012; Berrington et al., 2015). Both approaches have in common that they investigate whether and to which extent births can be ‘recuperated’ at later ages. In other words: Does the increase in the age at first birth simply shift the fertility schedule within the life course? Does it compress it into a shorter time period? Or does a late age at first childbearing eventually result in a decline in cohort fertility rates?

On the one hand, late age at first parenthood is assumed to result in lower completed fertility because of the rapid drop in fecundity with women's age (Billari et al., 2007). Albeit that the availability of assisted reproduction has fuelled a discussion on the biological limits of fertility, medical research usually points out that women's ability to conceive and bear children declines over time due to the ‘depletion and ageing of the pool of ovocytes stored in both ovaries during the fetal period’ (Velde et al., 2012, p. 1179). The ability to conceive and bear children declines gradually, but this process is assumed to accelerate around age 35. For men, fecundity is not subject to the same mechanisms, but the ability to father a child is nevertheless assumed to decay over time due to the decline of sperm quality with age(Schmidt et al., 2012; Johnson et al., 2015).

On the other hand, a late onset of fertility may not necessarily lead to lower completed fertility, because couples anticipate their impaired possibilities to have children at later ages. Postponement may be irrelevant for completed fertility, because couples have the ability to influence the spacing of births and may, thus, ‘squeeze’ their children into a shorter time interval. This ‘time squeeze’ has been assumed to particularly explain the fertility behaviour of work-oriented women who postpone childbearing to advance in their career and accelerate childbearing at later ages (Kreyenfeld, 2002; Gerster et al., 2007; Bartus et al., 2013). Furthermore, it was pointed out that the correlation between age at first birth and completed fertility was not a solid evidence of a causal relationship. Women and men who postpone first birth might do so because they had originally planned to have fewer children. Thus, couples self-select themselves into late childbearing, depending on their fertility preferences. Furthermore, genetic factors, such as health impairment, are important and often unobservable factors that affect the onset as well as completion of fertility (Kohler et al., 1999; Rodgers et al., 2008).

There is a significant body of empirical studies that has examined the association of age at first birth and subsequent fertility. However, the studies radically differ in terms of data and method employed. (Schmidt et al. 2012) correlated the period TFR and the mean age at first birth for several European countries for 2007, but did not find a strong association between the two measures. (Velde et al. 2012) examined the impact of fertility postponement on involuntary childlessness and total fertility in six European countries. Based on a simulation with macro-level data, they show that total fertility would have been between 0.03 and 0.05 higher in 2008 if there had been no postponement since the 1980s. (Andersson et al. 2009) rely on Scandinavian register data to show that early onset of childbearing leads to a higher number of children among women born between the 1930s and 1950. The mean number of children is well above 2 children for women who were age 25 or younger when they had their first child. For women who postponed first birth beyond age 35, total fertility drops below 1.5 children per woman. The effect of age at first birth on total fertility is linear and strong. But it is also shown that patterns vary across birth cohorts and countries.

This article adds to the discussion on the association between first birth and completed fertility by employing the promotion time model (Yakovlev and Tsodikov, 1996; Tsodikov, 1998; Chen et al., 1999) initially motivated to analyse cancer data. It belongs to the class of cure survival models extending classical event history models by acknowledging that an unknown proportion of the studied population will never experience the event of interest. (Bremhorst et al. 2016) motivated their use in fertility when studying the effect of the educational levels of a woman and of her partner on second and third parity progression. The motivation for the promotion in a fertility context works as follows. Assume that at the onset of the process (i.e., directly after the previous birth), a woman has $N \sim P (θ)$ (Poisson distributed) possible decisive arguments to opt for an additional child. Let F be the proper c.d.f. of the (independent and identically distributed) latent times $Y_{1}, \dots, Y_{n}$ necessary for any of these n arguments to initiate a new pregnancy. Then, one can show that the population survival and density functions of the event time are

S_{p} (t | θ, F) = exp (- θ F (t)); f_{p} (t | θ, F) = θ f (t) S_{p} (t | θ, F),

(1.1)

where $f (t) = \frac{\partial F (t)}{\partial t}$ is the latent density. In particular, the probability of never becoming pregnant again (or ‘cured’ in the jargon of the survival literature) is

P [N = 0] = exp (- θ) = lim_{t \to + \infty} S_{p} (t | θ, F) .

(1.2)

The expressions of the survival and hazard functions of the susceptible population are, respectively, given by

S_{u} (t | θ, F) = \frac{S_{p} (t | θ, F) - exp (- θ)}{1 - exp (- θ)}; h_{u} (t | θ, F) = \frac{f_{p} (t | θ, F)}{S_{p} (t | θ, F) - exp (- θ)} .

(1.3)

Independent baseline covariates, denoted by $x$ (including a constant ‘1’ to multiply the intercept) and $z$ (without such a constant), may enter, for example, the model through a log-link on parameter $θ$ and through a Cox model for $F (t)$ , respectively, yielding

\begin{matrix} θ (x) = exp (η_{θ} (x)) ​ ​ ​ & with & ​ ​ η_{θ} (x) = α^{T} x \end{matrix}

(1.4)

\begin{matrix} F (t | z) = 1 - S_{0} (t)^{exp (η_{F} (z))} ​ ​ ​ & with & ​ ​ η_{F} (z) = β^{T} z . \end{matrix}

(1.5)

Although the Cox model specified in (1.5) describes the distribution of the time necessary for each possible latent argument to become decisive to initiate a new pregnancy, the main interest for a demographer is in the conditional effect of a given covariate (for fixed values of the others) on the ‘global timing’ of an additional birth. This ‘global’ interpretation can be achieved by conditioning on the number of latent decisive arguments. Indeed, assume for a moment that $x$ and $z$ are restricted to be the same binary characteristics, such as to be unemployed ( $x = z = 0$ ) or to be employed ( $x = z = 1$ ). Assume also that the number of possible decisive arguments for a new pregnancy at the beginning of the follow-up are known and happen to be equal for two given women, one employed and one unemployed. Then, a negative $β$ regression coefficient in (1.5) means that the employed woman will tend to wait longer to have an additional child. On the other hand, it should be noted that an inspection of the simultaneous credible region of the susceptible survival $S_{u} (t | θ, F)$ or hazard $h_{u} (t | θ, F)$ function in (1.3) is required to conclude on the unconditional (on $N$ ) effect of a given covariate on the ‘global timing’ of an additional birth for susceptible mothers.

Considering the promotion time model defined in (1.1) combined with regression models (1.4) and (1.5), Tsodikov (2002) proposed, in a frequentist framework, a nonparametric estimation of the baseline survival function $S_{0} (t)$ . Within the Bayesian paradigm, Yin and Ibrahim (2005) assumed a piece-wise exponential distribution, while (Bremhorst and Lambert 2016) opted for a flexible specification of $S_{0} (t)$ using P-splines (Eilers and Marx, 1996, 2010). For recent articles using or extending the promotion time model, we refer the interested reader to Liu and Shen (2009), Kim et al. (2009), Lopes and Bolfarine (2012) and Li and Lee (2017).

Since women's fecundity decays rapidly after age 35, we assume that age at first birth had a non-linear effect on fertility timing and quantum. Therefore, in this work, instead of assuming a linear effect of the continuous covariates in the regression parts $η_{θ} (\cdot)$ and $η_{F} (\cdot)$ in (1.4) and (1.5), double additive models based on Bayesian P-splines (Lang and Brezger, 2004) will be specified. Flexible modelling of continuous covariates were already considered in many different contexts, see (Hastie and Tibshirani 1990) for an early reference and, to cite a few, (Wood 2006, 2011) for an efficient software implementation and Lambert (2013) for their double use (to model jointly the conditional location and dispersion) in a semi-parametric Bayesian framework.

The article is organized as follows: Sections 2 and 3 define the flexible additive models using Bayesian P-splines for the continuous covariates and for the logarithm of the baseline hazard function, respectively. Bayesian inference techniques, including the specification of the prior distributions and the description of the Metropolis-within-Gibbs algorithm to sample from the joint posterior are described in Section 4. A simulation study assessing the accuracy of the new methodology is presented in Section 5, while Section 6 focuses on the analysis of the German Socio-Economic Panel (GSOEP) datasets. A discussion concludes the article.

2 Flexible regression model for covariate effects

Let $x^{b} = (x_{1}^{b}, \dots, x_{P}^{b})$ and $z^{b} = (z_{1}^{b}, \dots, z_{Q}^{b})$ be the binary covariate vectors influencing, respectively, the probability of having the event and its timing for susceptible subjects. Denoted by $x^{c} = (x_{1}^{c}, \dots, x_{R}^{c})$ and $z^{c} = (z_{1}^{c}, \dots, z_{S}^{c})$ , the sets of continuous covariates assumed to take values in $(- 1, 1)$ and having an impact on the probability of being cured and on the timing of the event for susceptible individuals, respectively. The assumption made on the support of the continuous covariates is not restrictive since any continuous covariate can be transformed to meet that requirement.

The regression models defined in (1.4) and (1.5) assume that the continuous covariates contribute linearly to the functions $η_{θ} (\cdot)$ and $η_{F} (\cdot)$ involved in the conditional probability of being cured and in the conditional distribution of the event time for susceptible subjects. These linearity assumptions can be relaxed by specifying double additive models for the effects of continuous covariates:

\begin{matrix} η_{θ} (x^{b}, x^{c}) & = & α_{0} + α^{T} x^{b} + \sum_{r = 1}^{R} g_{r}^{θ} (x_{r}^{c}), \end{matrix}

(2.1)

\begin{matrix} η_{F} (z^{b}, z^{c}) & = & β^{T} z^{b} + \sum_{s = 1}^{S} g_{s}^{F} (z_{s}^{c}) . \end{matrix}

(2.2)

The functions $g_{r}^{θ} (.)$ and $g_{s}^{F} (.)$ are specified using a linear combination of a large number (say, $L$ ) of cubic B-splines:

g_{r}^{θ} (x_{r}^{c}) = \sum_{l = 1}^{L} ϕ_{rl}^{θ} {\tilde{b}}_{l} (x_{r}^{c}); g_{s}^{F} (z_{s}^{c}) = \sum_{l = 1}^{L} ϕ_{sl}^{F} {\tilde{b}}_{l} (z_{s}^{c}),

(2.3)

where ${{\tilde{b}}_{1} (.), \dots, {\tilde{b}}_{L} (.)}$ denotes a cubic B-splines basis common to all covariates and associated to a predefined number of equidistant knots on $[- 1, 1]$ . To avoid identification issues with the intercept $α_{0}$ (cf. Eq. 2.1) and with the baseline distribution of the Cox model (cf. Eq. 2.2), each spline coefficients vector $(ϕ_{1}^{θ}, \dots, ϕ_{R}^{θ})$ and $(ϕ_{1}^{F}, \dots, ϕ_{S}^{F})$ is constrained to have its components summing to 0.

To ensure smoothness and a linear behaviour in the limiting case, the likelihood is combined with a second order roughness penalty on finite differences of adjacent B-spline parameters (Eilers and Marx, 1996, 2010). For example, a second-order penalty associated to the $r th$ continuous covariate having an impact on the probability of the event is written as $τ_{r}^{θ} \sum_{l} {(Δ^{2} ϕ_{rl}^{θ})}^{2} = τ_{r}^{θ} {(ϕ_{r}^{θ})}^{T} {(D_{2})}^{T} D_{2} ϕ_{r}^{θ}$ , where $τ_{r}^{θ}$ is the penalty parameter and $D_{2}$ is the difference matrix defined as follows:

D_{2} = [\begin{matrix} 1 & - 2 & 1 & 0 & . . . & 0 \\ 0 & 1 & - 2 & 1 & . . . & 0 \\ ⋮ & ⋮ & ⋱ & ⋱ & ⋱ & ⋮ \\ 0 & 0 & . . . & 1 & - 2 & 1 \end{matrix}]

3 Flexible specification of the baseline distribution

As suggested by Bremhorst and Lambert (2016) in the context of the promotion time model, the baseline distribution $S_{0} (t)$ (of the Cox model defined in (1.5)) can be specified in a flexible way through the logarithm of the baseline hazard function $log (h_{0} (t))$ :

log (h_{0} (t)) = \sum_{k = 1}^{K} ϕ_{k} b_{k} (t),

(3.1)

where ${b_{1} (.), \dots, b_{K} (.)}$ denotes a cubic B-splines basis associated to a predefined number of equidistant knots on $[0, t_{max}]$ , where $t_{max}$ is the upper bound of the follow-up interval.

Remember that the latent cumulative distribution function $F (t)$ in (1.1) is assumed to be proper. When a non(semi)’parametric model is proposed for the latent distribution, that assumption can be forced using the zero tail constraint (Taylor, 1995; Zeng et al., 2006) forcing $F (t_{max}) = 1$ and, thus, $S_{u} (t_{max} | θ, F) = 0$ , see Eq. (1.3). Accordingly, Bremhorst and Lambert (2016) suggest to fix the last spline coefficient to a large enough value such that $S_{0} (t)$ smoothly decreases to $0$ when approaching the maximum follow-up time, thereby translating that a susceptible subject must have experienced the event of interest by that time. However, when the population survival function reaches its minimum value long before $t_{max}$ (as can be revealed from the wide plateau in the right tail of the Kaplan–Meier estimate of the population survival, when the maximum censoring time is much larger than the maximum observed failure time), the estimation of the spline parameters can be numerically instable since no information is available to estimate the spline parameters supported by knots located in that part of the follow-up interval. Therefore, in this article, we suggest to set all the spline parameters associated to knots located over the maximum observed failure time (i.e., in the plateau of the Kaplan–Meier estimate) to an arbitrary large value (such as 10).

As for the flexible additive models defined in (2.3), a roughness penalty on finite difference of adjacent B-spline coefficients is used to force smoothness. However, for the baseline distribution, a third-order penalty is used to force a limiting quadratic behaviour of the logarithm of the baseline hazard for large penalty parameter values.

4 Bayesian inference

4.1 Likelihood

Let $𝕕_{i} = (t_{i}, δ_{i}, x_{i}^{b}, x_{i}^{c}, z_{i}^{b}, z_{i}^{c})$ be the set of the observable data on the $i th$ individual under study, where $t_{i}$ and $δ_{i}$ denote the observed survival time and the event indicator (i.e., $δ_{i} = 0$ , if the subject is right censored). Since the data are potentially right censored, the log-likelihood is given by

l (Φ | 𝕕) = \sum_{i = 1}^{I} δ_{i} log (h_{p} (t_{i} | Φ, 𝕕_{i})) + log (S_{p} (t_{i} | Φ, 𝕕_{i})),

(4.1)

where $Φ$ is the set of the model parameters.

4.2 Bayesian model

The frequentist roughness penalty introduced by Eilers and Marx (1996, 2010) is translated in a Bayesian framework by Lang and Brezger (2004) into a multivariate normal prior distribution for the spline parameters:

\begin{matrix} ϕ | τ & \sim & N_{K} (0, {(τ P_{3})}^{- 1}); \end{matrix}

(4.2)

\begin{matrix} ϕ_{r}^{θ} | τ_{r}^{θ} & \sim & N_{L} (0, {(τ_{r}^{θ} P_{2})}^{- 1}) \forall r = 1, \dots, R; \end{matrix}

(4.3)

\begin{matrix} ϕ_{s}^{F} | τ_{s}^{F} & \sim & N_{L} (0, {(τ_{s}^{F} P_{2})}^{- 1}) \forall s = 1, \dots, S, \end{matrix}

(4.4)

where $P_{d} = D_{d}^{T} D_{d} + ε I$ is a full rank matrix for some small quantity $ε$ ( $10^{- 6}$ , say). Lang and Brezger (2004) suggest a Gamma distribution $𝕘 (a, b)$ with mean $\frac{a}{b}$ as a prior for the penalty parameters. However, it can be shown that the shape of the estimated curve could be significantly influenced by the choice of the values of the hyperparameters $a$ and $b$ . Therefore, in this work and as in Bremhorst and Lambert (2016) and in Bremhorst et al. (2016), we suggest to use the robust roughness penalty prior distribution proposed by Jullion and Lambert (2007):

\begin{matrix} τ | λ \sim 𝕘 (\frac{ν}{2}, \frac{ν λ}{2}) & ; & λ \sim 𝕘 (a, b); \end{matrix}

(4.5)

\begin{matrix} τ_{r}^{θ} | λ_{r}^{θ} \sim 𝕘 (\frac{ν}{2}, \frac{ν λ_{r}^{θ}}{2}) & ; & λ_{r}^{θ} \sim 𝕘 (a, b) \forall r = 1, \dots, R; \end{matrix}

(4.6)

\begin{matrix} τ_{s}^{F} | λ_{s}^{F} \sim 𝕘 (\frac{ν}{2}, \frac{ν λ_{s}^{F}}{2}) & ; & λ_{s}^{F} \sim 𝕘 (a, b) \forall s = 1, \dots, S . \end{matrix}

(4.7)

Bremhorst and Lambert (2016) showed that if the value of the hyperparameters $a$ and $b$ are small enough (say, $a = b = 10^{- 4}$ ), the posterior distribution of $ν$ is close to a uniform and, therefore, the value of $ν$ has no noticeable impact on the shape of the estimated curve.

Independent normal distributions with a large variance $σ^{2}$ are used as priors for the regression parameters:

(α_{0}, α, β) \sim N_{1 + p + q} (0, σ^{2} I) .

(4.8)

4.3 Posterior sampling using MCMC

The logarithm of the joint posterior distribution is obtained using Bayes’ theorem, as the sum of (4.1), the log-likelihood and the logarithm of the prior density functions defined in (4.2)– (4.8)

It can be shown that only the conditional posterior distributions of the penalty parameters belong to a known family of distributions:

\begin{matrix} τ | ϕ, λ, 𝕕 \sim 𝕘 (\frac{ν + K}{2}, \frac{ν λ + ϕ^{T} P_{3} ϕ}{2}); \\ λ | τ, 𝕕 \sim 𝕘 (a + \frac{ν}{2}, b + \frac{ν τ}{2}); \\ τ_{r}^{θ} | ϕ_{r}^{θ}, λ_{r}^{θ}, 𝕕 \sim 𝕘 (\frac{ν + L}{2}, \frac{ν λ_{r}^{θ} + {(ϕ_{r}^{θ})}^{T} P_{2} ϕ_{r}^{θ}}{2}) & \forall & r = 1, \dots, R; \\ λ_{r}^{θ} | τ_{r}^{θ}, 𝕕 \sim 𝕘 (a + \frac{ν}{2}, b + \frac{ν τ_{r}^{θ}}{2}) & \forall & r = 1, \dots, R; \\ τ_{s}^{F} | ϕ_{s}^{F}, λ_{s}^{F}, 𝕕 \sim 𝕘 (\frac{ν + L}{2}, \frac{ν λ_{s}^{F} + {(ϕ_{s}^{F})}^{T} P_{2} ϕ_{s}^{F}}{2}) & \forall & s = 1, \dots, S; \\ λ_{s}^{F} | τ_{s}^{F}, 𝕕 \sim 𝕘 (a + \frac{ν}{2}, b + \frac{ν τ_{s}^{F}}{2}) & \forall & s = 1, \dots, S . \end{matrix}

Therefore, a Metropolis-within-Gibbs algorithm will be used to sample from the joint posterior. When subsets of the model parameters are correlated, as with the P-spline models defined in (2.3) and (3.1), reparametrizing the joint posterior distribution using an adequate estimation of the posterior correlation structure of the model parameters might improve the mixing of the posterior chains (Lambert, 2007). Such an estimation can be achieved using a non-linear optimizer, based, for example, on the Augmented Lagrangian method (see, e.g., Nocedal and Wright, 2006, Chapter 17), enabling to deal with linear constraints on the parameters.

The reparametrization can be practically computed in the multivariate Metropolis steps in the following way: let $\hat{Σ_{r}^{θ}}$ denote the estimated variance–covariance matrix of the spline coefficients vector $ϕ_{r}^{θ}$ (associated to the $r$ th continuous covariate influencing the probability of being cured) obtained using the Augmented Lagrangian method aforementioned and let $\hat{L_{r}^{θ}}$ denote its Cholesky decomposition. Then, compute $ϕ_{r}^{θ (prop)} = ϕ_{r}^{θ (m - 1)} + \hat{L_{r}^{θ}} z$ , where $ϕ_{r}^{θ (m - 1)}$ is the value of $ϕ_{r}^{θ}$ at iteration $(m - 1)$ of the Metropolis-within-Gibbs algorithm, $z \sim N_{L} (0, σ_{r}^{θ} I)$ and $σ_{r}^{θ}$ are tuned during the burn-in period to reach the target acceptance rate, 23% in a multivariate setting (Haario et al., 2001; Atchadè and Rosenthal, 2005). For identifiability purpose, recentre $ϕ_{r}^{θ (prop)}$ . Then, accept (resp. reject) $ϕ_{r}^{θ (prop)}$ and set $ϕ_{r}^{θ (m)} = ϕ_{r}^{θ (prop)}$ (resp. $ϕ_{r}^{θ (m)} = ϕ_{r}^{θ (m - 1)})$ using a Metropolis step.

To speed up the convergence of MCMC algorithm, the posterior mode of the joint posterior distribution, for fixed values of $λ, λ_{1}^{θ}, \dots, λ_{R}^{θ}, λ_{1}^{F}, \dots, λ_{S}^{F}$ , is used as initial values of the MCMC algorithm.

5 Simulation

Table 1:
Simulation results for $S = 500$ replicates: Coverage probabilities (in $%)$ of the 95% credible intervals for each regression parameter

$n = 500$ $n = 1 000$

Cure probability Timing Cure probability Timing

Cure (%) Censored (%) $α_{0}$ $α_{1}$ $α_{2}$ $β_{1}$ $β_{2}$ $α_{0}$ $α_{1}$ $α_{2}$ $β_{1}$ $β_{2}$

$23$ $28$ 98.4 94.0 94.6 93.2 93.2 99.2 95.2 94.8 94.6 94.2

$23$ $50$ 98.2 92.0 91.8 92.6 90.8 97.4 93.8 94.8 92.2 94.8

$48$ $53$ 98.2 93.0 94.0 93.2 92.2 99.4 95.2 96.2 93.8 94.6

$48$ $69$ 98.2 93.4 92.6 92.0 90.8 98.8 94.2 93.8 91.8 92.8

		$n = 500$	$n = 1 000$
$23$	$28$	98.4	94.0	94.6	93.2	93.2	99.2	95.2	94.8	94.6	94.2
$23$	$50$	98.2	92.0	91.8	92.6	90.8	97.4	93.8	94.8	92.2	94.8
$48$	$53$	98.2	93.0	94.0	93.2	92.2	99.4	95.2	96.2	93.8	94.6
$48$	$69$	98.2	93.4	92.6	92.0	90.8	98.8	94.2	93.8	91.8	92.8

The numerical performances of the proposed methodology are illustrated through a simulation study. For all the settings described hereafter, $S = 500$ replicates of sample size $n = 500$ and $n = 1 000$ were generated. The simulated datasets are analysed using the model presented in Sections 2 and 3 combined with the Bayesian inference techniques described in Section 4. The flexible additive models (cf. Eq. 2.3) are defined using a basis of $L = 10$ B-splines associated to equidistant knots on $(- 1, 1)$ (for the relocated and rescaled continuous covariates), while a basis with $K = 15$ B-splines was taken on the follow-up interval to model (the log hazard of) the baseline distribution (cf. Eq. 3.1). The results presented in Tables 1 and 2 and in Figures 1–6 are based on one MCMC chain of length $75 000$ , including a burn-in period of length $25 000$ . The $z$ -scores of the Geweke convergence diagnostics (not reported) were used to assess the convergence of the Markov chains (Geweke, 1992).

Two percentages of cured individuals were considered: $23 %$ and $48 %$ . The percentages of right-censored subjects among the susceptible sub-population are controlled by two censoring distributions: a Weibull distribution with mean 17.86 and standard deviation 6.49 and an Exponential distribution with mean 18.18. Both are truncated at $t_{max} = 25$ , the upper bound of the follow-up, leading to $5 %$ and more than $20 %$ of right-censored susceptible subjects, respectively. Since the baseline distribution is specified, in each setting, as a Weibull distribution with mean $8$ and standard deviation $4.18$ , the sufficiently long follow-up assumption is respected. Indeed, less than $0.1 %$ of the event occur after time $t_{max} = 25$ under this distribution.

Let $W_{1}, W_{2} \sim$ Bernoulli $(0.5)$ and $W_{3}, W_{4} \sim U (- 1, 1)$ be the four independent covariates having simultaneously an impact on the probability of the event and on its timing for the susceptible sub-population. For each setting, one defines $x^{b} = (W_{1}, W_{2}) = z^{b}$ , $(β_{1}, β_{2}) = (0.4, - 0.4)$ , $x^{c} = (W_{3}, W_{4}) = z^{c}$ , $f_{1}^{θ} (.) = g_{2}^{F} (.)$ and $f_{2}^{θ} (.) = g_{1}^{F} (.)$ , where $f_{1}^{θ} (x) = \frac{x}{2}$ and $f_{2}^{θ} (x) = \frac{sin (2 π x)}{1.5}$ . The values of $α_{0}$ and $α$ is tuned in each setting to reach the target percentage of cure subjects.

Figure 1 shows a negligible bias of the posterior medians (as estimators) for all regression parameters in each setting. The variability of the estimates increases with the percentage of right-censored subjects among the susceptible sub-population. Not surprisingly, an increase of the proportion of right-censored susceptible subjects is more influential on the precision of parameter estimates than a comparable increase in the percentage of cured units. The coverage probabilities of the $95 %$ credible intervals for all the regression parameters (see Table 1) are close to their nominal value except for the intercept $α_{0}$ in the regression model for the cure probability where the coverage probabilities are a bit too large.

Figures 2 and 5 illustrate the estimations of a linear effect of a continuous covariate on the probability of having the event and on its timing for the susceptible sub-population. In each setting and in both parts of the model, the linear trend is captured without bias by the flexible additive models. As expected, the accuracy of the estimates increases with the sample size and decreases when the percentage of right censoring increases.

The estimation of a sine-like effect of a continuous covariate on the probability of being cured is pictured in Figure 3. In each considered scenario, the additive model was able to recover the non-linear shape. However, when the percentage of right censoring is important, a small bias (becoming negligible when the sample size increases) appears in the area where the concavity of the curve changes, as illustrated on the fourth row of Figure 3). This issue is more pronounced when estimating the non-linear effect on the timing of the event for susceptible subjects, as shown in Figure 4. This can be explained by the decreasing information on the regression parameters in the Cox model for susceptible subjects when the percentage of right-censored units increases. Regarding accuracy, the same conclusions as for the estimation of the linear effects hold.

Table 2 reports the coverage probabilities of the 95% simultaneous credible intervals for the additive terms in the two regression models. The estimated coverages are all close to their nominal values except for the non-linear term $g_{1}^{F} (\cdot)$ involved in the regression model for the timing of the event for susceptible individuals, when information is very sparse (typically when the sample size and the percentage of events are simultaneously small), leading to an oversmoothed estimate. Table 2 also reveals that for replicates yielding a credible region that does not contain entirely the true additive term, the average proportion (percentage in brackets) of the ( $-$ 1,1) covariate support leading to a coverage of the true value is quite satisfactory.

Following Held (2004), simultaneous credible intervals for the unknown functionals are computed using MCMC. This method is particularly useful when the joint posterior distribution of the parameters of interest (here, the parameters in the additive terms involved in the regression parts) does not belong to a known multivariate family. Practically, such simultaneous credible intervals can be computed, for example, using the function credible.region in the R package bayesSurv (Garcĺa-Zattera et al., 2016).

Finally, the same conclusions as for the other function estimates hold for the estimation of the baseline distribution, pictured in Figure 6.

Table 2:

Simulation results for $S = 500$ replicates: Coverage probabilities (in $%$ ) of the 95% simultaneous credible regions. In brackets: average proportion of the ( $-$ 1,1) covariate support for which the true additive term belongs to the simultaneous credible region when the credible region does not contain entirely the true additive term. The true functions are defined as $g_{1}^{θ} (x_{1}^{c}) = 0.5 x_{1}^{c}$ ; $g_{1}^{F} (z_{1}^{c}) = sin (2 π z_{1}^{c})$ ; $g_{2}^{θ} (x_{2}^{c}) = sin (2 π x_{2}^{c})$ ); $g_{2}^{F} (z_{2}^{c}) = 0.5 z_{2}^{c}$

			Cure probability		Timing
n	Cure (%)	Censored (%)	$g_{1}^{θ} (x_{1}^{c})$	$g_{2}^{θ} (x_{2}^{c})$	$g_{1}^{F} (z_{1}^{c})$	$g_{2}^{F} (z_{2}^{c})$
4*500	$23$	$28$	96.6 (92.54)	96.0 (95.22)	91.4 (93.22)	96.0 (91.47)
	$23$	$50$	96.8 (91.11)	93.6 (93.25)	87.4 (92.30)	95.6 (88.60)
	$48$	$53$	97.8 (94.21)	93.8 (93.52)	84.6 (93.69)	96.6 (89.11)
	$48$	$69$	97.4 (93.57)	91.0 (92.77)	78.6 (91.23)	96.2 (90.97)
4* $1 000$	23	28	98.4 (89.43)	97.8 (94.12)	96.8 (96.61)	97.6 (93.20)
	$23$	$50$	96.6 (95.08)	96.8 (94.50)	93.4 (95.46)	96.6 (91.28)
	$48$	$53$	98.8 (94.36)	99.2 (97.89)	97.4 (94.99)	96.8 (94.12)
	$48$	$69$	98.2 (92.32)	95.6 (96.07)	90.4 (92.64)	96.8 (91.85)

Figure 1:

Simulation results for $S = 500$ replicates: Boxplots of the errors of the regression parameter estimates (the posterior medians) for each considered scenario. Grey (resp. white) boxes correspond to settings where a small (resp. high) percentage of the susceptible population are randomly right censored

Figure 2:

Simulation results: estimation of $g_{1}^{θ} (x_{1}) = \frac{x_{1}}{2}$ . $S = 500$ replicates (one grey curve per dataset) with sample size $n = 500$ (left) or $n = 1 000$ (right). Each row refers to a percentage of cured and right-censored individuals. The solid line corresponds to the true function and the dashed line is the point-wise median of the 500 estimated curves

Figure 3:

Simulation results: estimation of $g_{2}^{θ} (x_{2}) = \frac{sin (2 π x_{2})}{1.5}$ . $S$ = 500 replicates (one grey curve per dataset) with sample size $n = 500$ (left) or $n = 1 000$ (right). Each row refers to a percentage of cured and right-censored individuals. The solid line corresponds to the true function and the dashed line is the point-wise median of the 500 estimated curves

Figure 4:

Simulation results: estimation of $g_{1}^{F} (z_{1}) = \frac{sin (2 π z_{1})}{1.5}$ . $S = 500$ replicates (one grey curve per dataset) with sample size $n = 500$ (left) or $n = 1 000$ (right). Each row refers to a percentage of cured and right-censored individuals. The solid line corresponds to the true function and the dashed line is the point-wise median of the 500 estimated curves

Figure 5:

Simulation results: estimation of $g_{2}^{F} (z_{2}) = 0.5 z_{2}$ . $S = 500$ replicates (one grey curve per dataset) with sample size $n = 500$ (left) or $n = 1 000$ (right). Each row refers to a percentage of cured and right-censored individuals. The solid line corresponds to the true function and the dashed line is the point-wise median of the 500 estimated curves

Figure 6:

Simulation results: estimation of the baseline distribution $S_{0} (t)$ . $S = 500$ replicates (one grey curve per dataset) with sample size $n = 500$ (left) or $n = 1 000$ (right). Each row refers to a percentage of cured and right-censored individuals. The solid line corresponds to the true function and the dashed line is the point-wise median of the 500 estimated curves

6 Application

Table 3:
Descriptive statistics of the independent covariates for all mothers at risk of an additional birth

Continuous variables (in years) Second birth Third birth

3Age at 1st birth (mother) Mean 27.9 26.4

Median 28.0 26.0

std 4.8 4.7

Age at previous birth (partner) Mean 30.6 32.7

Median 30.0 32.0

std 5.5 5.1

Birth interval $1 st$ / $2 nd$ child Mean $-$ 3.6

Median $-$ 3.0

std $-$ 2.2

Education levels Frequency

3Mother Low 340 393

Medium 936 959

High 231 238

Partner Low 311 309

Medium 901 945

High 295 336

Continuous variables (in years)	Second birth	Third birth
3*Age at 1st birth (mother)	Mean	27.9	26.4
	Median	28.0	26.0
	std	4.8	4.7
Age at previous birth (partner)	Mean	30.6	32.7
Median	30.0	32.0
std	5.5	5.1
Birth interval $1 st$ / $2 nd$ child	Mean	$-$	3.6
Median	$-$	3.0
std	$-$	2.2
Education levels	Frequency
3*Mother	Low	340	393
	Medium	936	959
	High	231	238
Partner	Low	311	309
Medium	901	945
High	295	336

Data for this analysis comes from the GSOEP. The GSOEP is the longest standing household panel survey in Europe (Wagner et al., 2007). Interviews are conducted with all household members aged $17$ and older (and since recently also with younger respondents) and are repeated every year. The original sample included West German households and an oversample of foreigners. Since its initiation in $1984$ , the SOEP has been extended several times. One of the most significant extensions was the inclusion of an East German subsample in the year $1990$ when German reunification was ratified. For our analysis, we use data from wave $2014$ (Release $31.1$ ). We restrict the analysis to the period 1984–2014. Our focus is the analysis of second and third birth transitions. Thus, we omit childless women from the investigation. Furthermore, we restrict the analysis to women aged between 17–49 who became under risk of having a second or third child after they entered the panel study. By doing so, we disregard left-truncated observations. We have also decided to right censor the cases $15$ years after previous birth. Although our analysis mainly focuses on the mother, we also considered age and education of the co-residential partner. Omitted are all respondents who did not have a partner at previous birth. Furthermore, cases are censored at the end of a co-residential union. We also omitted all cases with missing information on education of the respondent or her partner. Finally, we restricted the sample to women who had their preceding child in Western Germany. This restriction is attributed to the large East–West differences in fertility behaviour that prevailed after reunification. Also migrant women were omitted, as migrant fertility would require a separate investigation. Finally, we dropped cases from the so-called high income sample so that the analyses are not distorted due to the oversampling of high-income groups. The final sample comprises $1 507$ one-child mothers and $1 590$ two-child mothers. Among the one-child mothers, $57.0 %$ had a second child by the end of the follow-up. Among the two-child mothers, $22.8 %$ had a third child.

The dependent variables are the progression to the second and third child with time measured since the preceding birth. Transition to second and third births are studied separately. For both studies, the Kaplan–Meier estimate of the survival function, illustrated on Figure 7, shows a plateau in the right tail of the distribution suggesting that an unknown proportion of the studied one-child (resp. two-child) women will never have a second (resp. third) child. The data are analysed using the extended promotion time model defined in (1.1) with (2.1) and (2.2) as expressions for the link functions $η_{θ} (.)$ and $η_{F} (.)$ combined with the Bayesian inference methodology presented in Section 4. The convergence of the resulting posterior chains (of length $150 000$ with a burn-in period of $50 000$ ) were assessed by the trace plots and using convergence diagnostics tools such as those proposed by Geweke (1992). Our key independent variable is the age of the mother at first birth. We furthermore control for partner's age at the previous birth. For third birth, we also consider the duration between the first and second births. We control for the education levels of the mother and of her partner described as being either low (less than a vocational degree), medium (a vocational training degree) or high (a university or college degree). Education is handled as a time-constant covariate and was fixed at its value at the occasion of the previous birth (see Table 3 for the descriptive statistics).

Figure 7:

Transition probabilities to second and third births (Kaplan–Meier survival estimates)

6.1 Second birth

In a first step, we investigate second birth behaviour. Our interest is mainly the effect of mother's age at first birth. The results from this covariate on the probability of having a second birth are presented on the first row of Figure 8. The left panel shows the estimated additive term $g^{θ} (age)$ in (2.1) quantifying the effect of mother's age at first birth on the probability of having a second child for given education levels and age of the partner at first birth, while the right panel illustrates the evolution with age at first birth of the probability of having a second child for a medium-educated mother having a medium-educated co-residential partner who had his first child at age 30. There is a clear negative non-linear effect: the older a woman at first birth, the lower the probability of having a second child with an accelerating decline when the first birth occurred after 30. This finding is compatible with medical studies that indicate a decline of fecundity over time and a rapid deterioration in the mid-30s. If that interpretation is correct, one would expect that this effect is less pronounced for males. The results for partner's age at previous (first) birth, pictured on the second row of Figure 8, might partially support this claim. However, although second birth probabilities still decline, in particular, if the partner was 40 or older at first birth, the sample size is too small to assert it with a $90 %$ confidence level.

Figure 9 depicts (conditionally on the number of latent decisive arguments to initiate a second pregnancy) the effect of age at first birth on the fertility timing of a susceptible woman (for given educations levels and age of the partner at first birth). When a woman decides to have a second child, the time to that event tends to get shorter for mothers around their mid-30s at first birth. This result matches the finding for fertility quantum. Having a first child around age 35 seems, thus, to be a critical age. Women who have postponed childbearing up to that age either tend to accelerate childbearing or to remain one-child mothers. Due to small sample sizes, confidence bounds are rather large so that we have to be cautious with the interpretation. This particularly pertains to the effect of male's age at first birth on fertility timing. We do not find much of an association here.

The results from the other covariates, namely, the education levels of the two partners, are displayed in Table 4. The results corroborate earlier investigations with the GSOEP that focussed on the effect of education on fertility in cure fraction models (Bremhorst et al., 2016). For given ages of the two partners, highly educated women display a significantly higher probability of having a second child compared to medium, or less-educated women. In addition, having a highly educated partner significantly increases that probability compared to other couples.

Table 4:
Transition to second birth: Posterior median and 95% (HPD) credible interval for the regression parameter associated to categorical covariates. Sample size: $n = 1 507$ one-child women of whom $859$ ( $57 %$ ) reported having had a second child before censoring

Probability Timing

Est ${HPD}_{95 %}$ Est ${HPD}_{95 %}$

Intercept $-$ 0.696 [ $-$ 1.478; $-$ 0.010] $-$ $-$

Mother's education (ref. Medium)

Low $-$ 0.147 [ $-$ 0.378; 0.069] 0.061 [ $-$ 0.231; 0.354]

High 0.367 [0.047; 0.644] $-$ 0.171 [ $-$ 0.561; 0.198]

Partner's education (ref. Medium)

Low 0.108 [ $-$ 0.126; 0.335] $-$ 0.233 [ $-$ 0.524; 0.066]

High 0.406 [0.146; 0.653] 0.020 [ $-$ 0.310; 0.336]

	Probability	Timing
Intercept	$-$ 0.696	[ $-$ 1.478; $-$ 0.010]	$-$	$-$
Mother's education (ref. Medium)
Low	$-$ 0.147	[ $-$ 0.378; 0.069]	0.061	[ $-$ 0.231; 0.354]
High	0.367	[0.047; 0.644]	$-$ 0.171	[ $-$ 0.561; 0.198]
Partner's education (ref. Medium)
Low	0.108	[ $-$ 0.126; 0.335]	$-$ 0.233	[ $-$ 0.524; 0.066]
High	0.406	[0.146; 0.653]	0.020	[ $-$ 0.310; 0.336]

Figure 8:

Transition to second birth: Estimation of the conditional effect of the mother's age at first birth (row 1) and of the conditional effect of the partner's age at first birth (row 2) on the probability to have a second child. Column 2: Estimation of the conditional probability of having a second child w.r.t. the mother's age at first birth (row 1) and w.r.t the partner's age at first birth (row 2) (for median values of the other covariates). Dark (resp. light) grey region corresponds to the $95 %$ (resp. $90 %$ ) simultaneous credible region. Sample size: $n = 1507$ one-child women of whom $859$ ( $57 %$ ) reported having had a second child before censoring

Figure 9:

Transition to second birth: Estimation of the conditional effect of the mother's age at first birth (left) and of the conditional effect of the partner's age at first birth (right) on the timing of a second birth for susceptible women. Dark (resp. light) grey region corresponds to the $95 %$ (resp. $90 %$ ) simultaneous credible region. Sample size: $n = 1507$ one-child women of whom $859$ ( $57 %$ ) reported having had a second child before censoring

6.2 Third birth

In a second step, we examined the effect of mother's age at first birth on third-birth behaviour (for a given time interval between the first two births, a given partner's age at second birth and given education levels of the partners). The result of the impact of mother's age at first birth on third-birth parity progression is illustrated on the first row of Figure 10. The left-hand graph shows a very strong conditional (linear) negative effect of age at first birth on the probability of having a third child. For median values of the continuous covariates (cf. Table 3), the right-hand graph highlighted that this probability is close to 0 for a medium-educated two-child woman having a medium-educated co-residential partner and having had her first child after 35. On the other hand, the top-left graph of Figure 11 surmises that the conditional (on the number of latent decisive arguments available at second birth to initiate a third pregnancy) effect of age at first birth on the timing of third birth for susceptible two-child women is positive, suggesting that women squeeze their further children into shorter birth intervals. In particular, women who had their first child after age 30 tend to accelerate third birth transition. However, a larger number of births will be needed to claim it with at a 90% confidence level. The second row of Figure 10 displays the conditional effect of the time elapsed between the first two births on the probability to have a third child for fixed values of the other independent variables (left) and the evolution with the time interval between the first two births of that probability for a medium-educated woman having had her first child at age 26 and living with a medium-educated partner having had his second child at age 32 (right). It suggests that the longer the time interval between the first two children, the smaller the probability of having a third child. Although no significant conditional effect was found for the partner's age at second birth, the estimated conditional effects on the probability and on the timing of having a third child are pictured on the third row of Figure 10 and on the bottom-left graph of Figure 11, respectively.

Table 5 reports the effect of mother's and partner's education levels. It corroborates the findings in Bremhorst et al. (2016): While the mother's education level does not significantly affect the third-birth probability, a U-shaped effect for the partner's education is suggested, with medium-educated partner, significantly less likely to progress to a third child.

Table 5:
Transition to third birth: Posterior median and 95% (HPD) credible interval for the regression parameter associated to categorical covariates. Sample size: $n = 1 590$ two-child women of whom $364$ ( $22.8 %$ ) reported having had a third child before censoring

Probability Timing

Est ${HPD}_{95 %}$ Est ${HPD}_{95 %}$

Intercept $-$ 2.197 [ $-$ 3.292; $-$ 1.173] $-$ $-$

Mother's education (ref. Medium)

Low 0.218 [ $-$ 0.060; 0.501] $-$ 0.279 [ $-$ 0.647; 0.079]

High 0.290 [ $-$ 0.131; 0.711] $-$ 0.082 [ $-$ 0.596; 0.435]

Partner's education (ref. Medium)

Low 0.435 [0.152; 0.728] $-$ 0.371 [ $-$ 0.754; $-$ 0.006]

High 0.704 [0.363; 1.038] 0.056 [ $-$ 0.373; 0.475]

	Probability	Timing
Intercept	$-$ 2.197	[ $-$ 3.292; $-$ 1.173]	$-$	$-$
Mother's education (ref. Medium)
Low	0.218	[ $-$ 0.060; 0.501]	$-$ 0.279	[ $-$ 0.647; 0.079]
High	0.290	[ $-$ 0.131; 0.711]	$-$ 0.082	[ $-$ 0.596; 0.435]
Partner's education (ref. Medium)
Low	0.435	[0.152; 0.728]	$-$ 0.371	[ $-$ 0.754; $-$ 0.006]
High	0.704	[0.363; 1.038]	0.056	[ $-$ 0.373; 0.475]

Figure 10:

Transition to third birth—Column 1: Estimation of the conditional effect of the mother's age at first birth (row 1); of the conditional effect of the time elasped between first and second births (row 2) and of the conditional effect of the partner's age at second birth (row 3) on the probability to have a third child. Column 2: Estimation of the conditional probability of having a third child w.r.t. the mother's age at first birth (row 1); w.r.t the time elapsed between first and second births (row 2) and w.r.t the partner's age at second birth (row 3) (for median values of the other covariates). Dark (resp. light) grey region corresponds to the $95 %$ (resp. $90 %$ ) simultaneous credible region. Sample size: $n = 1 590$ two-child women of whom $364$ ( $22.8 %$ ) reported having had a third child before censoring

Figure 11:

Transition to third birth: Estimation of the conditional effect of the mother's age at first birth (top-left); of the conditional effect the time interval between first and second births (top-right) and of the conditional effect of the partner's age at second birth (bottom-left) on the timing of a third birth for susceptible women. Dark (resp. light) grey region corresponds to the $95 %$ (resp. $90 %$ ) simultaneous credible region. Sample size: $n = 1 507$ one-child women of whom $859$ ( $57 %$ ) reported having had a second child before censoring

7 Discussion

This article has used cure survival models to estimate the effect of age at first parenthood on fertility progression in Germany. To reach that goal, we considered double additive models to specify in a flexible way the effect of continuous covariates in cure fraction models. With this tool, we tried to detect non-linearities in the effect of mother's age at first birth on the progression to a second and third child. We suspected that age at first birth had a non-linear effect on fertility timing and quantum, because women's fecundity decays rapidly after age 35. Those who start their ‘fertility career’ at that age need to accelerate childbearing. Otherwise, they will be unable to progress to another child. Our analysis supports this view. While a high age at first birth seems to accelerate parity progression (specially for third birth) for susceptible women, it reduces the overall chance of having further children. This important finding underscores the relevance of cure fraction models in fertility research. Standard event history models easily generate misleading results because they conflate the effect of timing and quantum (see Bremhorst et al., 2016, Section 4.3 for more details). Cure survival models overcome this shortcoming, but have unfortunately not diffused much in social science research.

An important limitation to our research objectives in Section 6 was the limited sample size. We had about $1 500$ women at risk of second birth and another $1 500$ at risk of third birth at our disposal. These moderate values might explain several of the reported non-significant results. It also justifies why we decided not to study effect heterogeneities. In particular, we were not in a position to detect significant interaction effects of education and age at first birth. This is unfortunate because earlier studies suggested that the ‘meaning’ of age at birth depends on the education level (Hoem, 1996; Kreyenfeld, 2002; Bartus et al., 2013). While highly educated women regularly postpone childbirth in order to advance in their career, women with a low level of education who postpone parenthood are a select group, possibly with impaired abilities to have children. Thus, one would assume that a high age at first birth has a stronger negative effect on the probability of having a second or third child for low-educated women than for highly educated ones.

In a future study based on register data, it could be interesting to investigate the presence of an interaction between education and age at first birth to evaluate the aforementioned hypothesis.

From a theoretical point of view, a possible extension would be to develop a statistical tool to test the linearity of the effects of continuous covariates. To reach that goal, non-linear perturbations could be added to the linear specification of the continuous covariates in (1.4) and (1.5). These non-linear perturbations would be described using P-splines with a first order penalty and a zero mean to force a flat limiting behaviour of the correction for large values of the corresponding penalty parameter. An inspection of the simultaneous credible region for the correction and of the posterior of the associated penalty parameter could then be used to assess whether a linear approximation of the target additive term makes sense.

Acknowledgments

The authors disclosed receipt of the following financial support for the research, authorship and/or publication of this article: The authors acknowledge financial support from IAP research network P7/06 of the Belgian Government (Belgian Science Policy) and from the contract Projet d'Actions de Recherche Concertées (ARC) 11/16-039 of the Communauté française de Belgique, granted by the Académie Universitaire Louvain.

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 disclosed receipt of the following financial support for the research, authorship and/or publication of this article: The authors acknowledge financial support from IAP research network P7/06 of the Belgian Government (Belgian Science Policy) and from the contract Projet dActions de Recherche Concertées (ARC) 11/16-039 of the Communauté française de Belgique, granted by the Académie Universitaire Louvain.

References

Andersson

Knudsen

Neyer

Teschner

Ronsen

Lappegard

SkredeKandVikat

(2009) Cohort fertility patterns in the Nordic countries. Demographic Research , 20, 313–352.

Atchade

Rosenthal

(2005) On adaptive Markov chain Monte Carlo algorithms. Bernoulli , 11, 815–828.

Bartus

Murinko

Szalma

Szel

(2013) The effect of education on second births in Hungary: A test of the time-squeeze, self-selection, and partner effect hypotheses. Demographic Research , 28, 1–32.

Berrington

Stone

Beaujouan

(2015) Educational differences in timing and quantum of childbearing in Britain: A study of cohorts born 1940–1969. Demographic Research , 33, 733–764.

Billari

Kohler

H-P

(2004) Patterns of low and lowest-low fertility in Europe. Population Studies , 8, 161–176.

Billari

Kohler

H-P

Andersson

Lundstrom

(2007) Approaching the limit: Long-term trends in late and very late fertility. Population and Development Review , 33, 149–170.

Bratti

Tatsiramos

(2012) The effect of delaying motherhood on the second childbirth in Europe. Journal of Population Economics , 25, 291–321.

Bremhorst

Kreyenfeld

Lambert

(2016) Fertility progression in Germany: An analysis using flexible non parametric cure survival models. Demographic Research , 35, 505–534.

Bremhorst

Lambert

(2016) Flexible estimation in cure survival models using Bayesian P-splines. Computational Statistics and Data Analysis , 93, 270–284.

10.

Chen

Ibrahim

Sinha

(1999) A new Bayesian model for survival data with a surviving fraction. Journal of the American Statistical Association , 94, 909–919.

11.

Dommaraju

(2009). Timing of first birth in India. Genus , 65, 81–101.

12.

Eilers

PHC

Marx

(1996) Flexible smoothing with B-splines and penalties (with discussion). Statistical Science , 11, 89–121.

13.

Eilers

PHC

Marx

(2010) Splines, knots, and penalties.Wiley Interdisciplinary Reviews: Computational Statistics, 2, 637–653.

14.

Frejka

Sardon

J-P

(2006) First birth trends in developed countries: Persisting parenthood postponement. Demographic Research , 15, 147–180.

15.

Garca-Zattera

Jara

Komarek

(2016). A flexible aft model for misclassified clustered interval-censored data. Biometrics , 72, 473–483.

16.

Gerster

Keiding

Knudsen

Strandberg-Larsen

(2007) Education and second birth rates in Denmark 1981–1994. Demographic Research , 17, 181–210.

17.

Geweke

(1992) Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In Bayesian Statistics 4 Edition , edited by Bernardo

, Berger

Dawid

Smith

AFM

, pages 169–193. Oxford: Clarendon Press.

18.

Goldstein

Sobotka

Jasilioniene

(2009) The end of ‘lowest-low’ fertility? Population and Development Review , 35, 663–699.

19.

Haario

Saksman

Tamminen

(2001) An adaptive Metropolis algorithm. Bernoulli , 7, 223–242.

20.

Hastie

Tibshirani

(1990) Generalized Additive Models . London: Chapman & Hall.

21.

Held

(2004) Simultaneous posterior probability statements from Monte Carlo output. Journal of Computational and Graphical Statistics , 13, 20–35.

22.

Hoem

(1996) The social meaning of the age at second birth for third-birth fertility: A methodological note on the need to sometimes respecify an intermediate cariable. Yearbook of Population Research in Finland , 33, 333–339.

23.

Johnson

Dunleavy

Gemmell

Nakagawa

(2015) Consistent agedependent declines in human semen quality: A systematic review and metaanalysis. Ageing Research Reviews , 19, 22–33.

24.

Jullion

Lambert

(2007) Robust specification of the roughness penalty prior distribution in spatially adaptive Bayesian P-splines models. Computational Statistics and Data Analysis , 51, 2542–2558.

25.

Kim

Chen

(2009) A new latent cure rate marker model for survival data. Annals of Applied Statistics , 3, 1124–1146.

26.

Kohler

H-P

Rodgers

Christensen

(1999) Is fertility behavior in our genes? Findings from a Danish twin study. Population and Development Review , 25, 253–288.

27.

Kreyenfeld

(2002) Time squeeze, partner effect or self-selection? An investigation into the positive effect of women’s education on second birth risks in West Germany. Demographic Research , 7, 15–48.

28.

Lambert

(2007) Archimedean copula estimation using Bayesian splines smoothing techniques. Computational Statistics and Data Analysis , 51, 6307–6330.

29.

Lambert

(2013) Nonparametric additive locationscale models for interval censored data. Statistics and Computing , 23, 75–90.

30.

Lang

Brezger

(2004) Bayesian P-splines. Journal of Computational and Graphical Statistics , 13, 183–212.

31.

Lesthaeghe

(1994) The second demographic transition in Western countries: An interpretation. In Gender and Family Change in Industrialized Countries, edited by K Oppenheim Mason and A-M Jensen, pages 17–62. Oxford: Clarendon Press.

32.

Lesthaeghe

(2010) The unfolding story of the second demographic transition. Population and Development Review , 36, 211–251.

33.

Lee

(2017) A latent promotion time cure rate model using dependent tail-free mixtures. Journal of the Royal Statistical Society, Serie A. Doi: 10.1111/rssa.12226

34.

Liu

Shen

(2009) A semi-parametric accelerate failure time cure model. Journal of the American Statistical Association , 104, 1168–1178.

35.

Lopes

Bolfarine

(2012) Random effects in promotion time cure rate models. Computational Statistics and Data Analysis , 56, 75–87.

36.

Marini

Hodsdon

(1981). Effects of the timing of marriage and first birth of the spacing of subsequent births. Demography , 18, 529–548.

37.

Nocedal

Wright

(Eds) (2006) Numerical Optimization. In Springer Series in Operations Research and Financial Engineering . New York, NY: Springer-Verlag.

38.

Rodgers

Kohler

McGue

Behrman

Petersen

Bingley

Christensen

(2008) Education and cognitive ability as direct, mediating, or spurious influences on female age at first birth: Behavior genetic models fit to Danish twin data. American Journal of Sociology , 114, S202–32.

39.

Schmidt

Sobotka

Bentzen

Nyboe Andersen

(2012) Demographic and medical consequences of the postponement of parenthood? Human Reproduction Update , 28, 29–40.

40.

Sobotka

(2004) Is lowest-low fertility in Europe explained by the postponement of childbearing? Population and Development Review , 30, 195–220.

41.

Sobotka

(2008) Overview Chapter 6: The diverse faces of the second demographic transition in Europe. Demographic Research , 19, 171–224.

42.

Sobotka

Zeman

Lesthaeghe

Frejka

(2011) Postponement and recuperation in cohort fertility: New analytical and projection methods and their application. European Demographic Research Papers, pages 1–86. URL https://pdfs.semanticsc holar.org/d66c/3ad5b1f4437abca3b0ad9b 91edd8987d7929.pdf (last accessed 6 June 2018).

43.

Taylor

(1995) Semi-parametric estimation in failure time mixture models. Biometrics , 51, 899–907.

44.

Tsodikov

(1998) A proportional hazard model taking account of long-term survivors. Biometrics , 54, 1508–1516.

45.

Tsodikov

(2002) Semi-parametric model of long- and short-term survival: An application to the analysis of breast cancer survival in Utah by age and stage. Statistic in Medicine , 21, 895–920.

46.

Velde

Habbema

Leridon

Eijkemans

(2012) The effect of postponement of first motherhood on permanent involuntary childlessness and total fertility rate in six European countries since the 1970s. Human Reproduction , 27, 1179–1183.

47.

Wagner

Frick

Schupp

(2007) The German Socio-Economic Panel Study (SOEP): Scope, evolution and enhancements. Schmollers Jahrbuch , 127, 139–1369.

48.

Wood

(2006) Generalized Additive Models: An Introduction with R . Boca Raton, FL: CRC Press.

49.

Wood

(2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 73, 3–36.

50.

Yakovlev

Tsodikov

(1996) Stochastic Models for Tumor of Latency and Their Biostatistical Applications . Singapore: World Scientific Publishing.

51.

Yin

Ibrahim

(2005) Cure rate models: A unified approach. Canadian Journal of Statistics , 33, 559–570.

52.

Zeng

Yin

Ibrahim

(2006) Semiparametric transformation models for survival data with a cure fraction. Journal of the American Statistical Association , 101, 670–684.

	Probability		Timing
	Est	${HPD}_{95 %}$	Est	${HPD}_{95 %}$
Intercept	$-$ 0.696	[ $-$ 1.478; $-$ 0.010]	$-$	$-$
Mother's education (ref. Medium)
Low	$-$ 0.147	[ $-$ 0.378; 0.069]	0.061	[ $-$ 0.231; 0.354]
High	0.367	[0.047; 0.644]	$-$ 0.171	[ $-$ 0.561; 0.198]
Partner's education (ref. Medium)
Low	0.108	[ $-$ 0.126; 0.335]	$-$ 0.233	[ $-$ 0.524; 0.066]
High	0.406	[0.146; 0.653]	0.020	[ $-$ 0.310; 0.336]