## Abstract

Understanding conception probabilities is important not only for helping couples to achieve pregnancy but also in identifying acute or chronic reproductive toxicants that affect the highly timed and interrelated processes underlying hormonal profiles, ovulation, libido, and conception during menstrual cycles. Currently, 2 statistical approaches are available for estimating conception probabilities depending upon the research question and extent of data collection during the menstrual cycle: a survival approach when interested in modeling time-to-pregnancy (TTP) in relation to women or couples' purported exposure(s), or a hierarchical Bayesian approach when one is interested in modeling day-specific conception probabilities during the estimated fertile window. We propose a biologically valid discrete survival model that unifies the above 2 approaches while relaxing some assumptions that may not be consistent with human reproduction or behavior. This approach combines both the survival and the hierarchical models allowing investigators to obtain the distribution of TTP and day-specific probabilities during the fertile window in a single model. Our model allows for the consideration of covariate effects at both the cycle and the daily level while accounting for daily variation in conception. We conduct extensive simulations and utilize the New York State Angler Prospective Pregnancy Cohort Study to illustrate our approach. We also provide the code to implement the model in R software in the supplemental section of the supplementary material available at Biostatistics online.

## INTRODUCTION

Fecundability is defined as the probability of recognized conception in a menstrual cycle among couples having regular unprotected intercourse (Gini, 1924). It is used as a measure of a couple's fecundity or biologic capacity for reproduction. Motivated by the needs of natural family planners, and the research and clinical communities desire to identify fecundity determinants or timing of medical intervention, 2 quantities of interest have emerged: the time-to-pregnancy (TTP) and the day-specific conception probabilities in a given menstrual cycle for a couple.

TTP is defined as the number of menstrual cycles it takes for a couple, having regular sexual intercourse without the use of contraception, to conceive. It has been used as a measure to study the effect of various exposures on fecundity, see, for instance, Baird and others (1986), Law and others (2005), and Buck Louis and others (2009). Weinberg and Gladen (1986) first proposed a beta-geometric model for the probability distribution of the TTP. Subsequently, Scheike and Jensen (1997) proposed a discrete survival model for the hazard of conception in a cycle j given cycle level covariates $Uij$, as

(1.1)
where $Ri$ denotes a subject-specific random effect. A particular feature of this model is that the hazard is related linearly to the covariates when transformed by a complementary log–log function. This is precisely the proportional hazards model with random effects for grouped data. The Scheike and Jensen (1997) model allows for the inclusion of cycle-varying covariates but cannot incorporate the effects of the day-level covariates, typically collected in prospective pregnancy studies. Due to the inherent biases and recall errors in retrospectively ascertained TTPs [see, for instance, Weinberg and others, 1994a; Cooney and others, 2009), the current emphasis is on prospectively designed studies where extensive information is collected from a couple on various time scales including the woman level (e.g., previous reproductive history), cycle level (e.g., biomarkers for stress), and daily level (e.g., intercourse behavior, various reproductive hormonal levels, smoking, caffeine and alcohol consumption, and other lifestyle factors). Another limitation of using (1.1) is that it does not account for “immaculate conception,” that is, the hazard for conception in a cycle should be zero if no intercourse occurs.

A question of considerable interest, in the presence of detailed daily level information, is the probability of conception due to a single act of intercourse on a given day relative to ovulation. This is known as the day-specific conception probability. This enables one to determine the “fertile window,” a quantity of considerable interest to couples planning or trying to avoid pregnancy. The days outside of the “fertile window” have a small chance of pregnancy. The estimation of day-specific conception probabilities is complicated by the fact that it is unusual for a sexually active couple to have a single act of intercourse during the fertile window. Additionally, only a proxy for ovulation is available as the gold standard requires direct visualization of the ovaries via laparoscopic or ultrasonographic techniques that are available only for women seeking medical or infertility treatment such as in vitro fertilization or intracytoplasmic sperm injection. Consequently, differing precision exists in such identification.

The original seminal work of Barrett and Marshall (1969) proposed

where $Xij=(Xij1,…,XijK)$ denotes the indicator vector of intercourse in the fertile window and $pk$ is the day-specific probability of conception. The 2 main assumptions underlying this model are as follows.

• (A1) Each day of intercourse has an independent effect (that is to say, each ejaculate of sperm acts independently) to fertilize the ovum.

• (A2) The daily probabilities of conception are constant across couples and cycles.

Since then, the literature has responded to the second assumption (A2) to a great extent starting with the work of Schwartz and others (1980) who noted that conception is not only dependent on the timing of the intercourse but also on several biological factors such as the penetrability of the mucus, the capacity of the ovum to be fertilized and the receptivity of the uterine lining for implantation. This led to the model

where ω is the probability of a cycle being viable. Modifications of this model include the Weinberg and others (1994c) model that accommodates cycle-specific covariates; the Zhou and Weinberg (1996) model that incorporates day-specific covariates; the Zhou and others (1996) model that accounts for the within-woman dependency; the Dunson and Zhou (2000) model that incorporates both within-woman dependency and a sterile fraction. Zhou and Weinberg (1996) proposed an expectation–maximization (EM) algorithm for the maximum likelihood estimation of parameters and used a sandwich estimator of the variance to adjust for within-woman dependency. However, their estimates of day-specific probabilities are biased downward, and the sandwich estimator is not valid because less fertile women contribute more cycles to the data. This class of models culminated in the work of Dunson and Stanford (2005), who note the weak identifiability of the cycle viability term ω. Consequently, they proposed a hierarchical model
(1.2)
where $ξi$ denotes a couple-specific random effect. Note that in this model, the effects of covariates $Zij=(Zij1,…,ZijK)$ (possibly day-level) on the probability of conception in a menstrual cycle is mediated only through the daily level probabilities of conception.

An assumption that is common to all day-specific models is (A1), namely that the ejaculates of sperm introduced by the intercourse acts on different days compete with each other independently in an attempt to fertilize the ovum. In particular, the independence assumption necessitates the relation

(1.3)
where $Ak$ denotes the event that a sperm from kth intercourse act fertilizes the ovum, and $Xk$ denotes the indicator of intercourse occurring on the kth day of the cycle. This may not be a reasonable biological assumption. Each intercourse act in the fertile window introduces a fresh ejaculate of sperm in the reproductive tract that can potentially fertilize the ovum. It is well known in the clinical literature (van Duijn and Freund, 1971; Tyler and others, 1985; Levin and others, 1986; Carlsen and others, 2004) that tremendous inter- and intravariability exists in sperm quality ranging from azoospermic (absence of sperm) to high-quality sperm. In fact, studies have shown that there is a 29.2% reduction in sperm concentration when ejaculatory frequency went from one to two episodes during a 7-day period prior to semen collection following 2 days of abstinence. For 3 or more ejaculations, the level was reduced by about 41% as compared to the one ejaculate group (Carlsen and others, 2004). Consequently, this implies that the number of sperm available for fertilization does not just vary from day to day (assuming one ejaculate per day) but also depends on how many intercourse acts occurred previously (that is, on the previous $Xk$'s). Sperm concentration and number per ejaculate may affect the probability of conception. These data suggest that (1.3) may not be reasonable. Although we focused on sperm concentration for the purpose of illustration, a host of other semen characteristics, such as motility and morphology also impact conception probability and are dependent upon previous intercourse pattern (Carlsen and others, 2004), supporting the need to relax the independence assumption.

Royston and Ferreira (1999) proposed an alternative approach that does not require the independence assumption in that they assume that in cycles with multiple intercourse acts only the most fertile one contributes to the probability of conception. Although this may be a reasonable approximation in some cases, sperm introduced on less optimal days in the fertile window can also compete to fertilize the ovum and should not be ignored, as normal appearing sperm may survive up to 5 or 6 days in the female reproductive tract.

In this paper, we propose a model for the hazard for conception in a menstrual cycle that takes into account the daily level intercourse acts appropriate for conception as well as the effects of covariates (daily, cycle, or couple level). This enables us to assess the effect of covariates directly on the survival (hazard) function for TTP. This model retains the ease of interpretation of the effects of covariates as in the discrete Cox model, while still allowing us to assess the more subtle daily level probabilities of conception in a cycle. Moreover, in our approach, we do not require the assumption of independence of acts of intercourse in fertilizing the ovum within the same cycle inherent in the day-specific conception models. Under the assumption (A1), Dunson (2003) has discussed assessing hazard based on day-specific models that predate the Dunson and Stanford (2005) model. In Section 2, we present our model and discuss its relation with existing models; in Section 3, we provide extensive simulations; and in Section 4, we show the application of our model to a prospective pregnancy cohort study with preconception enrollment of women who were followed through 12 menstrual cycles at risk for pregnancy.

## MODEL FOR CONCEPTION

We begin by introducing some notation. Let $Ti$ denote the TTP for couple $i,i=1,2,…,n$. As is usual in many time to event studies, $Ti$ is subject to right censoring ($τi$) and one observes ${Ti∧τi,δi=I(Ti≤τi)}$, where $I(·)$ denotes the indicator function. Let $Xij=(Xij1,Xij2,…,XijK)$ denote the intercourse indicators in the fertile window of jth cycle for the ith couple. Denote by $Uij,$ the cycle-level covariates. Further, denote by $λi(j)$, the hazard rate for the TTP of the ith couple.

Let $Cj$ be the event that the ovum is fertilized in the jth cycle, and $Ak$ be the event that a sperm from the kth intercourse fertilizes the ovum. Note that the fertilization of an ovum normally requires a sperm originating from one of the potential intercourse acts that the couple may have had in the fertile window. Then $Cj=∪k=1KAk$ (disjoint union of events). So,

Under the independence assumption (A1) of the events $Ak$,
To avoid requiring the independence assumption (A1), we mimic the mixing of ejaculates of sperm from different intercourse acts in the reproductive tract of a woman by using an arbitrary linear combination of intercourse acts in the fertile window. In other words, we weigh separately the intercourse acts on different days so as to discriminate between an intercourse act occurring on day k with that occurring on day $k′$ in the fertile window. These weights are estimated based on the observed sample. Furthermore, we propose to directly model conception in cycle j, given that conception has not occurred so far, by
Observe that $P(Cj|C1c,…,Cj−1c,Xij,Uij)$ is the hazard for conception in cycle j. In other words, we propose the following discrete survival model for TTP:
(2.4)
We assume that the random effects, $exp(Ri)$, follow a Gamma distribution with mean 1 and variance η. Observe that the proposed model corrects for “immaculate conception,” that is, the hazard for conception in a cycle is zero if the couple does not have any intercourse in the fertile window of that cycle. The regression coefficients $αk$ capture the baseline kth day effect of intercourse on the probability of conception in cycle j. The cycle-varying parameter $ρ(j)$ denotes the cycle-specific baseline, a quantity of considerable interest (Weinberg and others, 1994b). The regression coefficients β capture the effect of the covariates $Uij$. Observe that if a couple had intercourse only on day k, that is, $Xijk=1,Xijk′=0,fork′≠k$, then, under the proposed model, the probability of conception in cycle j is given by
(2.5)
This is the probability of conception in cycle j if the couple had intercourse only on a specific day in the fertile window of cycle j. This is analogous to the day-specific probabilities of conception in the day-specific models for conception. Consequently, we refer to it as the kth day-specific conditional hazard of conception. Also, note that the effect of covariates on $λi(k)$ can be viewed as additive effects on complementary log–log scale.

Furthermore, we can also estimate the effect of covariates directly on the probability for conception in cycle j as follows:

Consequently, the probability mass function can be expressed as
This yields the survival function for $Ti$ as follows:
(2.6)
Consequently, using the proposed model (2.1), one can model the day-level hazard for conception via (2.2) as well as model the effects of covariates on the survival function (2.3) in the same model. Also interesting to note is the constant ratio of log survival functions for the time-independent covariates (assuming equal time-dependent covariates).

Similar to Scheike and Jensen (1997), the marginal forms (with respect to random effect $Ri$) for the probability of conception in cycle j and the hazard for conception in cycle j can be expressed as

Under the assumption of a Gamma distribution for $exp(Ri)$ with mean 1 and variance η, the hazard rate $λi(j)=P(Ti=j|Ti≥j)$ is given by
Thus, the marginal kth day-specific conditional probability of conception in cycle j is given by

Most prospective pregnancy studies design data collection to include the use of daily diaries to ascertain daily level covariates such as menstruation and sexual intercourse and, possibly, factors purported to impact couple fecundity (e.g., cigarette smoking, alcohol, and caffeine consumption). One can also estimate the effect of such covariates in the model by viewing them as a vector $Zij=(Zij1,…,ZijL)$, the number of days L of interest need not be the same as the fertile window. One can incorporate the daily-level covariate into the model as follows:

(2.7)
The observed likelihood for the discrete survival model (2.1) is similar to that of binary data with probability of success $λi(j).$ However, this is different from the Dunson and Stanford (2005) approach, where they pose the problem in terms of binary outcomes and model the covariate effects only through day-specific probabilities (1.2) coupled with the independence assumption (A1). Additionally, one cannot incorporate the cycle-varying intercept $ρ(j)$ in (1.2) due to lack of identifiability with day-varying baseline intercept, $αk.$ As mentioned previously, the effect $ρ(j)$ may be important in these prospective pregnancy studies.

To summarize, our model unifies the 2 approaches for modeling fecundability: TTP approach and day-specific approach while (i) accounting for the couple-level heterogeneity through the random effect, (ii) accounting for the cycle-level baseline effect $ρ(j)$, (iii) assessing day-level covariates directly on the cycle-level probability of conception rather than through day-specific probabilities of conception, (iv) while not requiring the independence of sperm fertilizing assumption (A1). The proposed discrete survival model can be fitted using a likelihood-based approach. We include the code to implement it using R software in the Supplementary Section of the supplementary material available at Biostatistics online.

## NUMERICAL STUDY

The goal of this section is to investigate the performance of the estimates using a likelihood-based approach. We also investigated the effect of zero-risk sets on estimation. Here, “zero-risk” indicates that a couple did not have an intercourse acts during the fertile window of a cycle and, consequently, did not put themselves at risk for conception. Observe that this unique requirement of an intermediate event is one of the features that sets the TTP data apart from the classical discrete survival setup. A common practice of incorporating intercourse into the discrete survival setup is to summarize the number of intercourse acts in the fertile window and include it as a cycle-varying covariate. Obviously, this practice ignores the effect of zero-risk sets on hazard. So, the simulation study focuses on the performance of the proposed estimates and also studies the impact of ignoring zero-risk sets.

We generated the data as follows: for subject i, the frailty variable $exp(Ri)$ was generated from a gamma distribution with mean 1 and variance η, the day-level intercourse behavior were generated such that $Xijk=0/1$ with probability of one as $0.6$ for all $k=1,…,6$, $j=1,…,12$, and the covariates $Ui=(Ui1,Ui2)$ were generated such that $Ui1=0/1$ variable with probability of one as $0.5$ each and $Ui2∼N(0,1)$. Subjects who had not experienced an event at $j=12$ were censored. Note that in this setup, we have incorporated approximately 8% zero-risk sets which is close to the percentage encountered in the real data analysis in the next section.

Tables 1 and 2 present the performance of the likelihood-based estimators accounting for the zero-risk sets as introduced in (2.1) or ignoring its effect. We made the comparison at various sample sizes ranging from $n=100,200,500$ to even $n=1000$ motivated by some recently completed prospective pregnancy studies, for example, the Life Study, the Oxford Conception Study. The censoring percentage was also varied from 10% to 30%. The results presented are based on 1000 replicates. The column bias refers to the average of the difference between the estimated value of the parameter and the true value, the Avg(SE) refers to the average of the asymptotic standard deviation (calculated using the estimated observed fisher information), SE(est) refers to the standard deviation of the estimates and CP refers to the coverage probability for a 95% confidence interval.

Table 1.

Summary of the simulation study with 10% censoring, $β1=−1,β2=0.5$, $α=(0.1,0.2,0.2,0.3,$$0.4,0.1)$, and $η=0.5$.

 Variable Corrected for zero-risk periods Not corrected for zero-risk periods Bias Avg(SE) SE(est) CP Bias Avg(SE) SE(est) CP $n=100$ $β1$ −0.01 0.30 0.32 0.942 0.00 0.30 0.32 0.937 $β2$ 0.01 0.16 0.16 0.944 0.01 0.16 0.16 0.938 $α1$ 0.00 0.25 0.25 0.950 0.17 0.25 0.25 0.884 $α2$ 0.00 0.25 0.26 0.952 0.17 0.25 0.26 0.877 $α3$ 0.00 0.25 0.26 0.940 0.17 0.25 0.26 0.892 $α4$ 0.01 0.25 0.25 0.954 0.18 0.25 0.25 0.889 $α5$ 0.01 0.25 0.26 0.955 0.19 0.25 0.26 0.867 $α6$ −0.01 0.25 0.25 0.949 0.16 0.25 0.25 0.898 η −0.02 0.21 0.23 0.910 −0.03 0.21 0.23 0.909 $n=200$ $β1$ −0.01 0.21 0.22 0.950 0.01 0.21 0.21 0.954 $β2$ 0.01 0.11 0.11 0.946 0.00 0.11 0.11 0.950 $α1$ 0.00 0.18 0.18 0.945 0.17 0.17 0.18 0.820 $α2$ −0.01 0.17 0.18 0.948 0.16 0.17 0.18 0.855 $α3$ 0.00 0.17 0.18 0.956 0.17 0.17 0.18 0.817 $α4$ 0.00 0.17 0.18 0.937 0.17 0.17 0.18 0.833 $α5$ 0.01 0.17 0.18 0.949 0.19 0.17 0.18 0.800 $α6$ 0.00 0.18 0.17 0.955 0.16 0.17 0.18 0.842 η −0.01 0.15 0.15 0.931 −0.02 0.15 0.15 0.935 $n=500$ $β1$ −0.01 0.13 0.14 0.946 0.01 0.13 0.14 0.948 $β2$ 0.00 0.07 0.07 0.958 0.00 0.07 0.07 0.958 $α1$ 0.00 0.11 0.11 0.961 0.16 0.11 0.11 0.699 $α2$ 0.00 0.11 0.11 0.939 0.17 0.11 0.11 0.661 $α3$ 0.00 0.11 0.11 0.944 0.16 0.11 0.11 0.679 $α4$ 0.00 0.11 0.11 0.944 0.17 0.11 0.11 0.639 $α5$ −0.01 0.11 0.11 0.948 0.17 0.11 0.11 0.666 $α6$ −0.01 0.11 0.11 0.945 0.15 0.11 0.11 0.696 η 0.00 0.09 0.10 0.948 −0.01 0.09 0.10 0.940 $n=1000$ $β1$ −0.01 0.09 0.09 0.949 0.01 0.09 0.09 0.948 $β2$ 0.00 0.05 0.05 0.967 −0.01 0.05 0.05 0.959 $α1$ 0.00 0.08 0.08 0.936 0.16 0.08 0.08 0.456 $α2$ 0.00 0.08 0.07 0.955 0.16 0.08 0.07 0.411 $α3$ 0.00 0.08 0.08 0.951 0.16 0.08 0.08 0.438 $α4$ 0.00 0.08 0.08 0.950 0.16 0.08 0.08 0.420 $α5$ 0.01 0.08 0.07 0.959 0.18 0.08 0.07 0.342 $α6$ 0.00 0.08 0.08 0.953 0.16 0.08 0.08 0.451 η 0.00 0.07 0.07 0.950 0.00 0.07 0.07 0.948
 Variable Corrected for zero-risk periods Not corrected for zero-risk periods Bias Avg(SE) SE(est) CP Bias Avg(SE) SE(est) CP $n=100$ $β1$ −0.01 0.30 0.32 0.942 0.00 0.30 0.32 0.937 $β2$ 0.01 0.16 0.16 0.944 0.01 0.16 0.16 0.938 $α1$ 0.00 0.25 0.25 0.950 0.17 0.25 0.25 0.884 $α2$ 0.00 0.25 0.26 0.952 0.17 0.25 0.26 0.877 $α3$ 0.00 0.25 0.26 0.940 0.17 0.25 0.26 0.892 $α4$ 0.01 0.25 0.25 0.954 0.18 0.25 0.25 0.889 $α5$ 0.01 0.25 0.26 0.955 0.19 0.25 0.26 0.867 $α6$ −0.01 0.25 0.25 0.949 0.16 0.25 0.25 0.898 η −0.02 0.21 0.23 0.910 −0.03 0.21 0.23 0.909 $n=200$ $β1$ −0.01 0.21 0.22 0.950 0.01 0.21 0.21 0.954 $β2$ 0.01 0.11 0.11 0.946 0.00 0.11 0.11 0.950 $α1$ 0.00 0.18 0.18 0.945 0.17 0.17 0.18 0.820 $α2$ −0.01 0.17 0.18 0.948 0.16 0.17 0.18 0.855 $α3$ 0.00 0.17 0.18 0.956 0.17 0.17 0.18 0.817 $α4$ 0.00 0.17 0.18 0.937 0.17 0.17 0.18 0.833 $α5$ 0.01 0.17 0.18 0.949 0.19 0.17 0.18 0.800 $α6$ 0.00 0.18 0.17 0.955 0.16 0.17 0.18 0.842 η −0.01 0.15 0.15 0.931 −0.02 0.15 0.15 0.935 $n=500$ $β1$ −0.01 0.13 0.14 0.946 0.01 0.13 0.14 0.948 $β2$ 0.00 0.07 0.07 0.958 0.00 0.07 0.07 0.958 $α1$ 0.00 0.11 0.11 0.961 0.16 0.11 0.11 0.699 $α2$ 0.00 0.11 0.11 0.939 0.17 0.11 0.11 0.661 $α3$ 0.00 0.11 0.11 0.944 0.16 0.11 0.11 0.679 $α4$ 0.00 0.11 0.11 0.944 0.17 0.11 0.11 0.639 $α5$ −0.01 0.11 0.11 0.948 0.17 0.11 0.11 0.666 $α6$ −0.01 0.11 0.11 0.945 0.15 0.11 0.11 0.696 η 0.00 0.09 0.10 0.948 −0.01 0.09 0.10 0.940 $n=1000$ $β1$ −0.01 0.09 0.09 0.949 0.01 0.09 0.09 0.948 $β2$ 0.00 0.05 0.05 0.967 −0.01 0.05 0.05 0.959 $α1$ 0.00 0.08 0.08 0.936 0.16 0.08 0.08 0.456 $α2$ 0.00 0.08 0.07 0.955 0.16 0.08 0.07 0.411 $α3$ 0.00 0.08 0.08 0.951 0.16 0.08 0.08 0.438 $α4$ 0.00 0.08 0.08 0.950 0.16 0.08 0.08 0.420 $α5$ 0.01 0.08 0.07 0.959 0.18 0.08 0.07 0.342 $α6$ 0.00 0.08 0.08 0.953 0.16 0.08 0.08 0.451 η 0.00 0.07 0.07 0.950 0.00 0.07 0.07 0.948

Avg(SE), the average of the asymptotic standard deviation; SE(est), the standard deviation of the estimates and CP, the coverage probability for a 95% confidence interval.

Table 2.

Summary of the simulation study with 30% censoring, $β=−1$, $α=(0.1,0.2,0.2,0.3,$$0.4,0.1)$, and $η=0.5$

 Corrected for zero-risk periods Not corrected for zero-risk periods Variable Bias Avg(SE) SE(est) CP Bias Avg(SE) SE(est) CP $n=100$ $β1$ −0.01 0.32 0.33 0.947 -0.01 0.32 0.33 0.949 $β2$ 0.01 0.16 0.17 0.944 0.01 0.16 0.17 0.946 $α1$ −0.01 0.27 0.27 0.953 0.12 0.27 0.27 0.927 $α2$ −0.01 0.27 0.28 0.954 0.12 0.26 0.28 0.913 $α3$ 0.00 0.27 0.27 0.950 0.14 0.26 0.26 0.927 $α4$ 0.01 0.26 0.27 0.947 0.15 0.26 0.27 0.905 $α5$ 0.01 0.26 0.25 0.963 0.15 0.26 0.25 0.915 $α6$ −0.01 0.27 0.27 0.957 0.12 0.27 0.27 0.918 η −0.02 0.31 0.34 0.900 −0.01 0.32 0.34 0.899 $n=200$ $β1$ 0.00 0.22 0.23 0.947 0.00 0.22 0.23 0.944 $β2$ 0.01 0.11 0.12 0.933 0.01 0.11 0.12 0.937 $α1$ 0.00 0.19 0.18 0.957 0.13 0.19 0.18 0.907 $α2$ 0.00 0.19 0.18 0.961 0.13 0.18 0.18 0.886 $α3$ 0.01 0.18 0.19 0.934 0.14 0.18 0.19 0.867 $α4$ 0.00 0.18 0.19 0.945 0.14 0.18 0.19 0.876 $α5$ 0.01 0.18 0.18 0.959 0.15 0.18 0.18 0.853 $α6$ 0.00 0.19 0.18 0.956 0.13 0.19 0.18 0.889 η −0.02 0.22 0.23 0.935 −0.02 0.22 0.23 0.934 $n=500$ $β1$ 0.00 0.14 0.15 0.935 0.00 0.14 0.15 0.935 $β2$ 0.00 0.07 0.07 0.959 0.00 0.07 0.07 0.955 $α1$ 0.00 0.12 0.12 0.958 0.13 0.12 0.12 0.786 $α2$ $−0.01$ 0.12 0.11 0.965 0.12 0.12 0.11 0.808 $α3$ $−0.01$ 0.12 0.12 0.949 0.13 0.12 0.12 0.797 $α4$ 0.01 0.12 0.12 0.948 0.15 0.11 0.12 0.735 $α5$ 0.00 0.11 0.12 0.948 0.14 0.11 0.12 0.751 $α6$ −0.01 0.12 0.12 0.948 0.12 0.12 0.12 0.795 η 0.00 0.14 0.14 0.954 0.00 0.14 0.14 0.952 $n=1000$ $β1$ 0.00 0.10 0.10 0.951 0.00 0.10 0.10 0.950 $β2$ 0.00 0.05 0.05 0.965 0.00 0.05 0.05 0.964 $α1$ 0.00 0.08 0.08 0.956 0.13 0.08 0.08 0.659 $α2$ 0.00 0.08 0.08 0.943 0.13 0.08 0.08 0.622 $α3$ 0.00 0.08 0.08 0.944 0.13 0.08 0.08 0.608 $α4$ 0.00 0.08 0.08 0.950 0.13 0.08 0.08 0.621 $α5$ 0.00 0.08 0.08 0.949 0.15 0.08 0.08 0.566 $α6$ 0.00 0.08 0.08 0.945 0.13 0.08 0.08 0.646 η 0.00 0.10 0.10 0.943 0.00 0.10 0.10 0.948
 Corrected for zero-risk periods Not corrected for zero-risk periods Variable Bias Avg(SE) SE(est) CP Bias Avg(SE) SE(est) CP $n=100$ $β1$ −0.01 0.32 0.33 0.947 -0.01 0.32 0.33 0.949 $β2$ 0.01 0.16 0.17 0.944 0.01 0.16 0.17 0.946 $α1$ −0.01 0.27 0.27 0.953 0.12 0.27 0.27 0.927 $α2$ −0.01 0.27 0.28 0.954 0.12 0.26 0.28 0.913 $α3$ 0.00 0.27 0.27 0.950 0.14 0.26 0.26 0.927 $α4$ 0.01 0.26 0.27 0.947 0.15 0.26 0.27 0.905 $α5$ 0.01 0.26 0.25 0.963 0.15 0.26 0.25 0.915 $α6$ −0.01 0.27 0.27 0.957 0.12 0.27 0.27 0.918 η −0.02 0.31 0.34 0.900 −0.01 0.32 0.34 0.899 $n=200$ $β1$ 0.00 0.22 0.23 0.947 0.00 0.22 0.23 0.944 $β2$ 0.01 0.11 0.12 0.933 0.01 0.11 0.12 0.937 $α1$ 0.00 0.19 0.18 0.957 0.13 0.19 0.18 0.907 $α2$ 0.00 0.19 0.18 0.961 0.13 0.18 0.18 0.886 $α3$ 0.01 0.18 0.19 0.934 0.14 0.18 0.19 0.867 $α4$ 0.00 0.18 0.19 0.945 0.14 0.18 0.19 0.876 $α5$ 0.01 0.18 0.18 0.959 0.15 0.18 0.18 0.853 $α6$ 0.00 0.19 0.18 0.956 0.13 0.19 0.18 0.889 η −0.02 0.22 0.23 0.935 −0.02 0.22 0.23 0.934 $n=500$ $β1$ 0.00 0.14 0.15 0.935 0.00 0.14 0.15 0.935 $β2$ 0.00 0.07 0.07 0.959 0.00 0.07 0.07 0.955 $α1$ 0.00 0.12 0.12 0.958 0.13 0.12 0.12 0.786 $α2$ $−0.01$ 0.12 0.11 0.965 0.12 0.12 0.11 0.808 $α3$ $−0.01$ 0.12 0.12 0.949 0.13 0.12 0.12 0.797 $α4$ 0.01 0.12 0.12 0.948 0.15 0.11 0.12 0.735 $α5$ 0.00 0.11 0.12 0.948 0.14 0.11 0.12 0.751 $α6$ −0.01 0.12 0.12 0.948 0.12 0.12 0.12 0.795 η 0.00 0.14 0.14 0.954 0.00 0.14 0.14 0.952 $n=1000$ $β1$ 0.00 0.10 0.10 0.951 0.00 0.10 0.10 0.950 $β2$ 0.00 0.05 0.05 0.965 0.00 0.05 0.05 0.964 $α1$ 0.00 0.08 0.08 0.956 0.13 0.08 0.08 0.659 $α2$ 0.00 0.08 0.08 0.943 0.13 0.08 0.08 0.622 $α3$ 0.00 0.08 0.08 0.944 0.13 0.08 0.08 0.608 $α4$ 0.00 0.08 0.08 0.950 0.13 0.08 0.08 0.621 $α5$ 0.00 0.08 0.08 0.949 0.15 0.08 0.08 0.566 $α6$ 0.00 0.08 0.08 0.945 0.13 0.08 0.08 0.646 η 0.00 0.10 0.10 0.943 0.00 0.10 0.10 0.948

Avg(SE), the average of the asymptotic standard deviation; SE(est), the standard deviation of the estimates and CP, the coverage probability for a 95% confidence interval.

Observe that the bias of the estimates of the regression coefficients for the zero-risk set corrected model are reasonably small even for small sample size, and their performance improves considerably as n increases. Furthermore, the estimated standard deviation gets closer to the sampling standard deviation with the sample size. The coverage probability for the estimates are close to the nominal 95% and become closer as the sample size increases. However, the considerable bias in estimates of the $αk$'s when the zero-risk sets are ignored causes the coverage probabilities to be much lower than the nominal value. In addition, the estimates of β and η have higher bias and standard error when zero-risk sets are ignored for both settings and all n. Overall, the estimates for the corrected model perform very well even for small sample size. We further observe that the coverage probability for the 95% confidence interval for the variance of the random effect was lower than that for other parameters for $n=100.$ However, it becomes much closer to the nominal level of 95% as the sample size increases from $n=200$ to $n=500$ and $n=1000.$

We next present our analysis of New York State Angler Prospective Pregnancy Cohort Study.

## ANALYSIS OF NEW YORK ANGLER PROSPECTIVE PREGNANCY COHORT STUDY

We illustrate our proposed method by analyzing the New York State Angler Prospective Pregnancy Cohort Study (Buck Louis and others, 2009). This prospective cohort study recruited women aged 20 to 34 years from 16 counties surrounding Lakes Erie and Ontario who were discontinuing contraception for the purposes of becoming pregnant. Women were followed until a human chorionic gonadotropin detected pregnancy was observed or up to 12 menstrual cycles at risk for pregnancy. A nice feature of this study is the follow-up time of 12 “at-risk” menstrual cycles, which is much longer than the 6 cycle follow-up of most other studies, see Buck and others (2004). Note that clinically a couple is eligible for infertility treatment if they do not conceive by 12 menstrual cycles. Among the 113 women recruited, 14 were pregnant at baseline and, thereby, excluded. Eighty-three (84%) women completed daily diaries on menstruation, sexual intercourse, home pregnancy test results and covariates believed to impact female fecundity. Given the absence of a proxy biomarker for ovulation in the study, we utilized the Ogino-Knaus method for estimating ovulation by counting back 14 days from the end of the cycle (Knaus, 1929; Ogino, 1930). A priori, the fertile window was defined as comprising eight days before ovulation through three days after ovulation (Buck Louis and others, 2009). In our analysis, we will focus on the following covariates, namely, female age (years) upon enrollment, parity (yes/no), cigarette smoking (yes/no) during cycle. We fitted the model

where $Uij=(Parityi,Agei,Smokeij).$ Here, we have assumed that $exp(Ri)$ follows a Gamma distribution with mean 1 and variance $η.$ The estimates, standard errors and 95% confidence intervals for the effect of Parity, Age, Smoking and the day-level intercourse behavior are given in Table 3.

Table 3.

Estimates from the proposed models of the New York State Angler Prospective Pregnancy Cohort Study

 Variable Estimate SE 95% CI Parity 1.235 0.424 (0.404, 2.065) Age 0.237 0.174 (−0.105, 0.578) Smoke −1.079 0.454 (−1.969, −0.189) $α−8$ −0.611 0.381 (−1.357, 0.136) $α−7$ 0.403 0.336 (−0.256, 1.062) $α−6$ 0.348 0.337 (−0.314, 1.009) $α−5$ 0.393 0.329 (−0.252, 1.039) $α−4$ 1.088 0.331 (0.439, 1.737) $α−3$ 0.511 0.300 (−0.077, 1.099) $α−2$ −0.285 0.323 (−0.918, 0.347) $α−1$ −0.375 0.333 (−1.027, 0.277) $α0$ −0.356 0.336 (−1.014, 0.302) $α+1$ −0.497 0.329 (−1.143, 0.148) $α+2$ 0.404 0.304 (−0.193, 1.000) $α+3$ −0.477 0.335 (−1.134, 0.181) η 0.132 0.374 (−0.601, 0.866)
 Variable Estimate SE 95% CI Parity 1.235 0.424 (0.404, 2.065) Age 0.237 0.174 (−0.105, 0.578) Smoke −1.079 0.454 (−1.969, −0.189) $α−8$ −0.611 0.381 (−1.357, 0.136) $α−7$ 0.403 0.336 (−0.256, 1.062) $α−6$ 0.348 0.337 (−0.314, 1.009) $α−5$ 0.393 0.329 (−0.252, 1.039) $α−4$ 1.088 0.331 (0.439, 1.737) $α−3$ 0.511 0.300 (−0.077, 1.099) $α−2$ −0.285 0.323 (−0.918, 0.347) $α−1$ −0.375 0.333 (−1.027, 0.277) $α0$ −0.356 0.336 (−1.014, 0.302) $α+1$ −0.497 0.329 (−1.143, 0.148) $α+2$ 0.404 0.304 (−0.193, 1.000) $α+3$ −0.477 0.335 (−1.134, 0.181) η 0.132 0.374 (−0.601, 0.866)

In Figure 1, we present the conditional day-specific probabilities of conception in cycle j for a non-smoking woman aged 30 years (the mean), and with parity = 1 versus 0. Note that the Figure 1(a) indicates the probability of conception in the first cycle given that the couple had intercourse on a specific day ($k=−8,⋯,+3$) in the fertile window. Figure 1(b) indicates the probability of conception in the second cycle given that the couple had intercourse only on a specific day in the fertile window and had not conceived in the first cycle, (c) and (d) denote these probabilities given that the couple had not conceived by the second and third cycles, respectively. These plots do indicate that parity significantly increases the conception probabilities at the daily-level. Also, the graphs do indicate considerable measurement error associated with identifying the day of ovulation by the Ogino-Knauss method, which is consistent with findings of using other proxies of ovulation with similar measurement errors. This is acknowledged in the reproductive literature, see Lynch and others (2006) for a summary.

Fig. 1.

Plot of $λ(k)(j)$ for the New York State Angler Prospective Pregnancy Cohort Study for an average-age nonsmoking couple by parity for (a) cycle 1, (b) cycle 2, (c) cycle 3, and (d) cycle 4.

Fig. 1.

Plot of $λ(k)(j)$ for the New York State Angler Prospective Pregnancy Cohort Study for an average-age nonsmoking couple by parity for (a) cycle 1, (b) cycle 2, (c) cycle 3, and (d) cycle 4.

Figure 2 displays the effect of parity on the TTP survival distribution for a nonsmoking woman aged 30 years with parity = 1 versus 0 and having intercourse on days −4 through 2. Here, it is clear that the parity significantly reduces the TTP, where the estimated median TTP reduces from 11 to 3 cycles.

Fig. 2.

Plot of $S^(t)$ for the New York State Angler Prospective Pregnancy Cohort Study by parity.

Fig. 2.

Plot of $S^(t)$ for the New York State Angler Prospective Pregnancy Cohort Study by parity.

Finally, we also assessed the model fit based upon a comparison with (i) model with no information concerning intercourse (ii) model with just accounting for total number of intercourse in the fertile window and ignoring the issue of “zero-risk” set as is the common practice (iii) proposed model incorporating the day-specific intercourse behavior in the fertile window. The Akaike information criterion for models (i–iii) were 355.17, 348.42, and 344.78, respectively. Thus, indicating improved model fit when accounting for daily intercourse pattern.

## DISCUSSION

We have proposed a unified approach for modeling fecundity, combining the discrete survival model for TTP and modeling the day-specific probabilities of conception while relaxing the sperm independence assumption. This unified approach is biologically plausible and provides an accessible method to ascertain the effects of covariates directly on the hazard for conception in a cycle. The day-specific conditional probabilities of conception have a practical interpretation for couples who have not yet conceived with regard to their chances in the next cycle, given their current pattern of intercourse behavior. In fact, this approach may be desirable when a prospective pregnancy study does not have daily level hormonal measurements either due to the burden on the subjects or the cost involved in obtaining such information. In such situations, the “fertile window” may be identified with considerable measurement error. Consequently, studying the effects of covariates through the day-level probabilities of conception approach of the hierarchical models may not be the best strategy. We also have illustrated the performance of the proposed likelihood-based estimates in the Angler study. Our findings are indicative of measurement error in specifying the day of ovulation, given the absence of a proxy biomarker in this cohort study. This limitation underscores the importance of cautious interpretation of the estimated day-specific probabilities of conception for this prospective cohort study. To this end, our findings await corroboration from future cohort studies such as the recently completed National Institute of Child Health and Human Development LIFE Study that utilized reliable proxies of ovulation for better-fit models.

An advantage of our approach is that it allows modeling TTP using conventional methods such as the proportional odds or discrete Cox models while accounting for day-level intercourse and other covariates to ensure biologically plausible models for estimation. Our approach also provides a context for assessing the sterile fraction in the context of the so-called “cure fraction.” Finally, another aspect of considerable interest is to assess how the biologically meaningful relaxation of assumption (A1) translates quantitatively. Note, the Dunson and Stanford (2005) model works under the independence assumption (A1), while our model does not have this requirement. An interesting problem to consider in future is to develop a model to account for various dependence structures and use the model that gives the best fit.

## SUPPLEMENTARY MATERIAL

Supplementary material is available at http://biostatistics.oxfordjournals.org.

## FUNDING

Intramural Research Program of the National Institutes of Health; Eunice Kennedy Shriver National Institute of Child Health and Human Development.

The authors would like to acknowledge that this study utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, MD (http://biowulf.nih.gov). Conflict of Interest: None declared.

## References

Baird
DDAY
Wilcox
AJ
Weinberg
CR
Use of time to pregnancy to study environmental exposures
American Journal of Epidemiology
,
1986
, vol.
124
pg.
470

Barrett
JC
Marshall
J
The risk of conception on different days of the menstrual cycle
Population Studies
,
1969
, vol.
23
(pg.
455
-
461
)
Buck
GM
Lynch
CD
Stanford
JB
Sweeney
AM
Schieve
LA
Rockett
JC
Selevan
SG
SM
Prospective pregnancy study designs for assessing reproductive and developmental toxicants
Environmental Health Perspectives
,
2004
, vol.
112
pg.
79

Buck Louis
GM
Dmochowski
J
Lynch
C
Kostyniak
P
McGuinness
BM
Vena
JE
Polychlorinated biphenyl serum concentrations, lifestyle and time-to-pregnancy
Human Reproduction
,
2009
, vol.
24
(pg.
451
-
458
)
Carlsen
E
Petersen
JH
AM
Skakkebaek
NE
Effects of ejaculatory frequency and season on variations in semen quality
Fertility and Sterility
,
2004
, vol.
82
(pg.
358
-
366
)
Cooney
MA
Buck Louis
GM
Sundaram
R
McGuiness
BM
Lynch
CD
Validity of self-reported time to pregnancy
Epidemiology
,
2009
, vol.
20
pg.
56

Dunson
DB
Incorporating heterogeneous intercourse records into time to pregnancy models
Mathematical Population Studies
,
2003
, vol.
10
(pg.
127
-
143
)
Dunson
DB
Stanford
JB
Bayesian inferences on predictors of conception probabilities
Biometrics
,
2005
, vol.
61
(pg.
126
-
133
)
Dunson
DB
Zhou
H
A Bayesian model for fecundability and sterility
Journal of the American Statistical Association
,
2000
, vol.
95
(pg.
1054
-
1062
)
Gini
C
Premierès recherces sur la fécondabilité de la femme
Proceedings of the International Mathematical Congress
,
1924
, vol.
2
(pg.
889
-
892
)
Knaus
H
Ovulationstermin und Konzeptionstermin
Zentralblatt Für Gynäkologie
,
1929
, vol.
53
(pg.
464
-
479
)
Law
DCG
Klebanoff
MA
Brock
JW
Dunson
DB
Longnecker
MP
Maternal serum levels of polychlorinated biphenyls and 1, 1-dichloro-2, 2-bis (p-chlorophenyl) ethylene (DDE) and time to pregnancy
American Journal of Epidemiology
,
2005
, vol.
162
pg.
523

Levin
RM
Latimore
J
Wein
AJ
Van Arsdalen
KN
Correlation of sperm count with frequency of ejaculation
Fertility and Sterility
,
1986
, vol.
45
(pg.
732
-
734
)
Lynch
CD
Jackson
LW
Buck Louis
GM
Estimation of the day-specific probabilities of conception: current state of the knowledge and the relevance for epidemiological research
Paediatric and Perinatal Epidemiology
,
2006
, vol.
20
(pg.
3
-
12
)
Ogino
K
Eine neue Methods zur Bestimmung des Ovulationstermines
Zentralblatt Für Gynäkologie
,
1930
, vol.
54
pg.
2193

Royston
P
Ferreira
A
A new approach to modeling daily probabilities of conception
Biometrics
,
1999
, vol.
55
(pg.
1005
-
1013
)
Scheike
TH
Jensen
TK
A discrete survival model with random effects: an application to time to pregnancy
Biometrics
,
1997
, vol.
53
(pg.
318
-
329
)
Schwartz
D
Macdonald
PDM
Heuchel
V
Fecundability, coital frequency and the viability of ova
Population Studies
,
1980
, vol.
34
(pg.
397
-
400
)
Tyler
JP
Matson
PL
Crockett
NG
The effects of frequent ejaculation on seminal spermatozoal number and calculation of daily spermatozoal output
Clinical Reproduction and Fertility
,
1985
, vol.
3
(pg.
145
-
149
)
van Duijn
C
Jr.
Freund
M
The relationship between some seminal characteristics and ejaculation frequency in the human male
European Journal of Obstetrics and Gynecology
,
1971
, vol.
5
(pg.
167
-
174
)
Weinberg
CR
Baird
DD
Wilcox
AJ
Sources of bias in studies of time to pregnancy
Statistics in Medicine
,
1994
, vol.
13

5–7
(pg.
671
-
681
)
Weinberg
CR
Florack
EIM
Zielhuis
GA
Rolland
R
The influence of occupational activity on the menstrual cycle and fecundability
Epidemiology
,
1994
, vol.
5
(pg.
476
-
477
)
Weinberg
CR
BC
The beta-geometric distribution applied to comparative fecundability studies
Biometrics
,
1986
, vol.
42
(pg.
547
-
560
)
Weinberg
CR
BC
Wilcox
AJ
Models relating the timing of intercourse to the probability of conception and the sex of the baby
Biometrics
,
1994
, vol.
50
(pg.
358
-
367
)
Zhou
H
Weinberg
CR
Modeling conception as an aggregated Bernoulli outcome with latent variables via the EM algorithm
Biometrics
,
1996
, vol.
52
(pg.
945
-
954
)
Zhou
H
Weinberg
CR
Wilcox
AJ
Baird
DD
A random effects model for cycle viability in fertility studies
Journal of the American Statistical Association
,
1996
, vol.
91
(pg.
1413
-
1422
)