PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of blackwellopenThis ArticleFor AuthorsLearn MoreSubmit
Statistics in Medicine
 
Stat Med. 2017 September 20; 36(21): 3361–3379.
Published online 2017 May 22. doi:  10.1002/sim.7345
PMCID: PMC5575519

The forecasting of menstruation based on a state‐space modeling of basal body temperature time series

Abstract

Women's basal body temperature (BBT) shows a periodic pattern that associates with menstrual cycle. Although this fact suggests a possibility that daily BBT time series can be useful for estimating the underlying phase state as well as for predicting the length of current menstrual cycle, little attention has been paid to model BBT time series. In this study, we propose a state‐space model that involves the menstrual phase as a latent state variable to explain the daily fluctuation of BBT and the menstruation cycle length. Conditional distributions of the phase are obtained by using sequential Bayesian filtering techniques. A predictive distribution of the next menstruation day can be derived based on this conditional distribution and the model, leading to a novel statistical framework that provides a sequentially updated prediction for upcoming menstruation day. We applied this framework to a real data set of women's BBT and menstruation days and compared prediction accuracy of the proposed method with that of previous methods, showing that the proposed method generally provides a better prediction. Because BBT can be obtained with relatively small cost and effort, the proposed method can be useful for women's health management. Potential extensions of this framework as the basis of modeling and predicting events that are associated with the menstrual cycles are discussed. © 2017 The Authors. Statistics in Medicine Published by John Wiley & Sons Ltd.

Keywords: basal body temperature, menstrual cycle length, periodic phenomena, sequential prediction, state‐space model

1. Introduction

The menstrual cycle is the periodic changes that occur in the female reproductive system that make pregnancy possible. Throughout the menstrual cycle, basal body temperature (BBT) also follows a periodic pattern. The menstrual cycle consists of two phases, the follicular phase followed by the luteal phase, with ovulation occurring at the transition between the two phases. During the follicular phase, BBT is relatively low with the nadir occurring within 1 to 2 days of a surge in luteinizing hormone that triggers ovulation. After the nadir, the cycle enters the luteal phase and BBT rises by 0.3°C to 0.5°C 1.

Considerable attention has been paid to the development of methods to predict the days of ovulation and the day of onset of menstruation; currently available methods include urinary and plasma hormone analyses, ultrasound monitoring of follicular growth, and monitoring of changes in the cervical mucus or BBT. Monitoring the change in BBT is straightforward because it requires neither expensive instruments nor medical expertise. However, of the currently available methods, BBT measurement is the least reliable because it has an inherently large day‐to‐day variability 1; therefore, statistical analyses are required to improve predictions of events associated with the menstrual cycle based on BBT.

Many statistical models of the menstrual cycle length have been proposed, with the majority explaining the marginal distribution of menstrual cycle length, which is characterized by a long right tail. To explain within‐individual heterogeneity, which is an important source of variation in menstrual cycle length, Harlow and Zeger 2 made the biological assumption that a single menstrual cycle is composed of a period of ‘waiting’ followed by the ovarian cycle. They classified menstrual cycles into standard, normally distributed cycles that contained no waiting time and into nonstandard cycles that contained waiting time. Guo et al. 3 extended this idea and proposed a mixture model consisting of a normal distribution and a shifted Weibull distribution to explain the right‐tailed distribution. By accommodating covariates, they also investigated the effect of an individual's age on the moments of the cycle length distribution. Lin et al. 4 described a linear mixed model that can explain heterogeneous variances of within‐woman cycle length. Huang et al. 5 attempted to model the changes in the mean and variance of menstrual cycle length that occur during the approach of menopause. By using a change‐point model, they found changes in the annual rate of change of the mean and variance of menstrual cycle length, which indicate the start of early and late menopausal transition. In contrast, Bortot et al. 6 focused on the dynamic aspect of menstrual cycle length over time. By using a state‐space modeling approach, they derived a predictive distribution of menstrual cycle length that was conditional on past time series. By integrating a fecundability model into their time series model, they also developed a framework that estimates the probability of conception that is conditional on the within‐cycle intercourse behavior. A related joint modeling of menstrual cycle length and fecundability has been attained more recently 7, 8.

Although daily fluctuations in BBT are associated with the events of the menstrual cycle, and therefore BBT time series analyses likely provide information regarding the length of the current menstrual cycle, there have been no previous studies modeling BBT data to predict menstrual cycle length. Thus, in the present study, we developed a statistical framework that provides a predictive distribution of menstrual cycle length (which, by extension, is a predictive distribution of the next day of onset of menstruation) that is sequentially updated with daily BBT data. We used a state‐space model that includes a latent phase state variable to explain daily fluctuations in BBT and to derive a predictive distribution of menstrual cycle length that is dependent on the current phase state.

This paper is organized as follows: in Section 2, we briefly describe the data required for the proposed method and how the test dataset was obtained. In Section 3, we provide the formulation of the proposed state‐space model for menstrual cycle and give an overview of the filtering algorithms that we use to estimate the conditional distributions of the latent phase variable given the model and dataset. In the same section, we also describe the predictive probability distribution for the next day of onset of menstruation that was derived from the proposed model and the filtering distribution for the menstrual phase. In Section 4, we apply the proposed framework to a real dataset and compare the accuracy of the point prediction of the next day of onset of menstruation among the proposed method and previous methods. In Section 5, we discuss the practical utility of the proposed method with respect to the management of women's health and examine the potential challenges and prospects for the proposed framework.

2. Data

We assumed that for a given female subject, a dataset containing a daily BBT time series and days of onset of menstruation was available. The BBT could be measured with any device (e.g., a conventional thermometer or a wearable sensor), although different model parameters may be adequate for different measurement devices. In the application of the proposed framework described in Section 4, we used real BBT time series and menstruation onset data that was collected via a website called Ran's story (QOL Corporation, Ueda, Japan), which is a website that allows registered users to upload their self‐reported daily BBT and days of menstruation onset to QOL Corporation's data servers. At the time of registration to use the service, all users of Ran's story agree to the use of their data for academic research. Although no data regarding the ethnic characteristics of the users were available, it is assumed that the majority of, if not all, the users were ethnically Japanese because Ran's story is provided only in the Japanese language.

3. Model description and inferences

3.1. State‐space model of the menstrual cycle

Here we develop a state‐space model for a time series of observed BBT, y t, and an indicator of the onset of menstruation, z t, obtained for a subject for days t=1,…,T. By z t=1, we denote that menstruation started on day t, whereas z t=0 indicates that day t is not the first day of menstruation. We denote the BBT time series and menstruation data obtained until time t as Y t=(y 1,…,y t) and Z t=(z 1,…,z t), respectively.

We considered the phase of the menstrual cycle, θtR(t=0,1,,T), to be a latent state variable. We let ϵ t as the daily advance of the phase and assume that it is a positive random variable that follows a gamma distribution with shape parameter α and rate parameter β, which leads to the following system model:

θt=θt1+ϵt
(1)

ϵtGamma(α,β).
(2)

Under this assumption, the conditional distribution of θ t, given θ t−1, is a gamma distribution with a probability density function:

p(θt|θt1)=Gamma(α,β)=βαΓ(α)(θtθt1)α1expβ(θtθt1).
(3)

It is assumed that the distribution for the observed BBT y t is conditional on the phase θ t. As periodic oscillation throughout each phase is expected for BBT, a finite trigonometric series is used to model the average BBT. Assuming a Gaussian observation error, we expressed the observation model for BBT as

yt=a+m=1M(bmcos2mπθt+cmsin2mπθt)+et
(4)

etNormal(0,σ2),
(5)

where M is the maximum order of the series. Conditional on θ t,y t then follows a normal distribution with a probability density function:

p(yt|θt)=Normalμ(θt),σ2=12πσ2expytμ(θt)22σ2,
(6)

where μ(θt)=a+m=1M(bmcos2mπθt+cmsin2mπθt). By this definition, μ(θ t) is periodic in terms of θ t with a period of 1.

For the onset of menstruation, we assume that menstruation starts when θ t ‘steps over’ the smallest following integer. This is represented as follows:

zt=0whenθt=θt1
(7)

=1whenθt>θt1,
(8)

where [left floor]x[right floor] is the floor function that returns the largest previous integer for x. Writing this deterministic allocation in a probabilistic manner, which is conditional on (θ t,θ t−1),z t follows a Bernoulli distribution:

p(zt|θt,θt1)=(1zt)I(θt=θt1)+ztI(θt>θt1),
(9)

where I(x) is the indicator function that returns 1 when x is true or 0 otherwise.

Let ξ=(α,β,σ,a,b1,,bM,c1,cM) be a vector of the parameters of this model. Given a time series of BBT, Y T=(y 1,…,y T), an indicator of menstruation, Z T=(z 1,…,z T), and a distribution specified for initial states, p(θ 1,θ 0), these parameters can be estimated by using the maximum likelihood method. The log‐likelihood of this model is expressed as

l(ξ;YT,ZT)=logp(y1,z1|ξ)+t=2Tlogp(yt,zt|Yt1,Zt1,ξ),
(10)

where

logp(y1,z1|ξ)=logp(y1|θ1)p(z1|θ1,θ0)p(θ1,θ0)dθ1dθ0,
(11)

and for t=2,…,T,

logp(yt,zt|Yt1,Zt1,ξ)=logp(yt|θt)p(zt|θt,θt1)p(θt,θt1|Yt1,Zt1)dθtdθt1,
(12)

which can be sequentially obtained by using the Bayesian filtering technique described later. Note that although p(y t|θ t)(t[gt-or-equal, slanted]1) and p(θ t,θ t−1|Y t−1,Z t−1)(t[gt-or-equal, slanted]2) depend on ξ, this dependence is not explicitly described for notational simplicity.

3.2. State estimation and calculation of log‐likelihood by using the sequential Bayesian filtering

Given the state‐space model of menstrual cycle described earlier and its parameters and data, the conditional distribution of an unobserved menstrual phase can be obtained by using recursive formulae for the state estimation problem, which are referred to as the Bayesian filtering and smoothing equations 9. We describe three versions of this sequential procedure, that is, prediction, filtering, and smoothing, for the state‐space model described earlier.

Let p(θ t,θ t−1|Y t,Z t) be the joint distribution for the phase at successive time points t and t−1, which is conditional on the observations obtained by time t. This conditional distribution accommodates all the data obtained by time t and is called the filtering distribution. Similarly, the joint distribution for the phase of successive time points t and t−1 conditional on the observations obtained by time t−1,p(θ t,θ t−1|Y t−1,Z t−1), is referred to as the one‐step‐ahead predictive distribution. For t=1,…,T, these distributions are obtained by sequentially applying the following recursive formulae:

p(θt,θt1|Yt1,Zt1)=p(θt|θt1)p(θt1|Yt1,Zt1)=p(θt|θt1)p(θt1,θt2|Yt1,Zt1)dθt2,
(13)

p(θt,θt1|Yt,Zt)=p(yt,zt|θt,θt1)p(θt,θt1|Yt1,Zt1)p(yt,zt|θt,θt1)p(θt,θt1|Yt1,Zt1)dθtdθt1=p(yt|θt)p(zt|θt,θt1)p(θt,θt1|Yt1,Zt1)p(yt|θt)p(zt|θt,θt1)p(θt,θt1|Yt1,Zt1)dθtdθt1,
(14)

where for t=1, we set p(θ 1,θ 0|Y 0,Z 0) as p(θ 1,θ 0), which is the specified initial distribution for the phase. Equations (13) and (14) are the prediction and filtering equations, respectively. Note that the denominator in Equation (14) is the likelihood for data at time t(Equations (11) and (12)). Hence, the log‐likelihood of the state‐space model is obtained through application of the Bayesian filtering procedure. Missing observations can be handled in the filtering equation (14) by marginalizing likelihood over the missing values.

The joint distribution for the phase, which is conditional on the entire set of observations, p(θ t,θ t−1|Y T,Z T), is referred to as the (fixed‐interval type) smoothed distribution. With the filtering and the one‐step‐ahead predictive distributions, the smoothed distribution is obtained by recursively using the following smoothing formula:

p(θt,θt1|YT,ZT)=p(θt,θt1|Yt,Zt)p(θt+1,θt|YT,ZT)p(θt+1|θt)p(θt+1,θt|Yt,Zt)dθt+1=p(θt,θt1|Yt,Zt)p(θt+1,θt|YT,ZT)dθt+1p(θt|Yt,Zt)=p(θt,θt1|Yt,Zt)p(θt+1,θt|YT,ZT)dθt+1p(θt,θt1|Yt,Zt)dθt1.
(15)

In general, these recursive formulae are not analytically tractable for non‐linear, non‐Gaussian, state‐space models. However, the conditional distributions, as well as the log‐likelihood of the state‐space model, can still be approximated by using filtering algorithms for general state‐space models. We use Kitagawa's non‐Gaussian filter 10, where the continuous state space is discretized into equally spaced grid points at which the probability density is evaluated. We describe the numerical procedure to obtain these conditional distributions by using the non‐Gaussian filter in Appendix A.1.

3.3. Predictive distribution for the day of onset of menstruation

A predictive distribution for the day of onset of menstruation is derived from the assumptions of the model and the filtering distribution of the state.

We denote the distribution function and the probability density function of the gamma distribution with shape parameter s and rate parameter r as G(·;s,r) and g(·;s,r), respectively. Let the accumulated advance of the phase be denoted by Δk(t), which is then calculated as Δk(t)=r=t+1t+kϵr,k=1,2,. We consider the conditional probability that the next menstruation has occurred before day k+t given the phase state θ t. We denote this conditional probability as F(k|θt)=PrΔk(t)θtθt, where [left ceiling]x[right ceiling] is the ceiling function that returns the smallest following integer for x. Therefore, F(k|θ t) represents the conditional distribution function for the onset of menstruation. Under the assumption of the state‐space model described earlier, F(k|θ t) is given as

F(k|θt)=θtθtg(x;kα,β)dx=1G(θtθt;kα,β).
(16)

The conditional probability function for the day of onset of menstruation, denoted as f(k|θ t), is then given as

f(k|θt)=F(k|θt)F(k1|θt)=1G(θtθt;kα,β)1Gθtθt;(k1)α,β=Gθtθt;(k1)α,βG(θtθt;kα,β),
(17)

where we set F(0|θ t)=0.

The marginal distribution for the day of onset of menstruation, denoted as h(k|Y t,Z t), can also be obtained with the marginal filtering distribution for the phase state, p(θ t|Y t,Z t). It is given as

h(k|Yt,Zt)=f(k|θt)p(θt|Yt,Zt)dθt.
(18)

Note that this distribution is conditional on the data obtained by time t, and therefore, it provides a predictive distribution for the day of onset of menstruation that accommodates all the available information. The point prediction for the day of onset of menstruation can also be obtained from h(k|Y t,Z t) and a natural choice for it would be the k that gives the highest probability, maxh(k|Yt,Zt).

In Appendix A, we describe the numerical procedure for the non‐Gaussian filtering that we used to obtain these predictive distributions.

4. Application

We used the BBT time series and menstruation onset data provided by 20 users of the Ran's story website (QOL Corporation). Data were collected through the course of 44 to 91 consecutive menstrual cycles (see Table 1 for a summary of the data). Data of each subject's first 29 consecutive menstrual cycles were used to fit the state‐space model described earlier and to obtain maximum likelihood estimates of the parameters. For the order of the trigonometric series to be determined, models with M=1,2,…,12(referring to models M1, M2, …, and M12, respectively) were fitted and then compared based on the Akaike information criterion (AIC) for each subject. The best AIC model and the remaining menstrual cycle data (i.e., the data that were not used for the parameter estimations) were then used to evaluate the accuracy of the prediction of the next day of onset of menstruation based on the root mean square error (RMSE) and the mean absolute error (MAE) for each subject. We compared accuracy of the sequential prediction, obtained by using the proposed method, with that of prediction based on two previous methods: (i) the conventional calendar calculation method, which predicts the upcoming menstruation day as the day after a fixed number of days from the onset of preceding menstruation, and (ii) the sequential predictive method proposed by Bortot et al. 6, which utilizes the time series of past cycle length to yield a model‐based prediction of the length of the current cycle. The technical details for this comparison are described in Appendix B.1.

Table 1

Summary of self‐reported menstrual cycle data obtained from 20 subjects.

Parameter estimates in the best AIC model for the 20 subjects are shown in Table 2. The selected order of the trigonometric series ranged from 2 to 12. The estimated relationship between expected BBT and menstrual phase, and the estimated probability density distribution for phase advancement for each subject are shown in Figure Figure1.1. Although the estimated temperature‐phase regression lines were ‘squiggly’ owing to the high order of the trigonometric series, they in general exhibited two distinct stages: the temperature tended to be lower in the first half of the cycle and higher in the second half, which is consistent with the well‐known periodicity of BBT through the menstrual cycle 1. Figure Figure22 shows examples of the estimated conditional distributions for menstrual phase and the associated predictive distributions for the day of onset of menstruation.

Figure 1

Model components for each subject. Each color represents a different subject. Lines show the best Akaike information criterion model associated with the maximum likelihood estimates for each subject. (A) Relationship between expected temperature and menstrual ...

Figure 2

Example estimated conditional distributions for menstrual phase (top panels) and the associated predictive distributions for the day of onset of menstruation (bottom panels). In the top panels, predictive and filtering distributions are shown as dashed ...

Table 2

Results of fitting the state‐space model by using the maximum likelihood method.

The RMSE of prediction of the next day of onset of menstruation provided by each method is shown in Figure Figure33 and Table 3. In the proposed method, the RMSE tended to decrease as the day approached the next day of onset of menstruation, suggesting that accumulation of the BBT time series data contributed to increasing the accuracy of the prediction. Except for subjects 9 and 20, the proposed method exhibited a better predictive performance than the calendar calculation and the method of Bortot et al. 6, for at least one of the time points of prediction we considered (i.e., the day preceding the next day of onset of menstruation) (Figure (Figure3).3). Compared with the best prediction provided by the calendar calculation method, the range, mean, and median of the rate of the maximum reduction in the RMSE for each subject in the sequential method were 0.066–1.481, 0.581, and 0.488, respectively.

Figure 3

Root mean square error (RMSE) of the prediction of the day of onset of the next menstruation. Solid lines indicate the RMSE of the sequential prediction (based on the best Akaike information criterion model for each subject), and horizontal dashed lines ...

Table 3

Root mean square error of prediction of the next day of onset of menstruation based on the calendar calculation method.

The predictive performance tended to be even better when the accuracy was measured by using the MAE; the MAE of the predictions of the next day of onset of menstruation provided by each method is shown in Figure Figure44 and Table 4. Compared with the best prediction provided by the calendar calculation method, the range, mean, and median of the rate of the maximum reduction in the MAE for each subject by using the sequential method was 0.056–1.368, 0.461, and 0.361, respectively.

Figure 4

Mean absolute error (MAE) of the prediction of the day of onset of the next menstruation. Solid lines indicate the MAE of the sequential prediction (based on the best Akaike information criterion model for each subject), and horizontal dashed lines indicate ...

Table 4

Mean absolute error of prediction of the next day of onset of menstruation based on the calendar calculation method.

5. Discussion

Here we constructed a statistical framework that provides a model‐based prediction of the day of onset of menstruation based on a state‐space model and an associated Bayesian filtering algorithm. The model describes the daily fluctuation of BBT and the history of menstruation. The filtering algorithms yielded a filtering distribution of menstrual phase that was conditional on all of the data available at that point in time, which was used to derive a predictive distribution for the next day of onset of menstruation.

The predictive framework we developed has several notable characteristics that make it superior to previous methods for predicting the day of onset of menstruation. State‐space modeling and Bayesian filtering techniques enable the proposed method to yield sequential predictions of the next day of onset of menstruation based on daily BBT data. Even though menstrual cycle length fluctuates stochastically, the prediction of day of onset of menstruation is automatically adjusted based on the daily updated filtering distribution, yielding a flexible yet robust prediction of the next day of onset of menstruation. This is markedly different from the conventional calendar calculation method that yields only a fixed prediction that is never adjusted. In a study similar to the present study, Bortot et al. 6 constructed a predictive framework for menstrual cycle length by using the state‐space modeling approach. However, in their model, the one‐step‐ahead predictive distribution for menstrual cycle length is obtained based on the past time series of cycle length, and within‐cycle information is not taken into account in the prediction.

Furthermore, compared with previous methods, the proposed method generally yielded a more accurate prediction of day of onset of menstruation. As the day approaches the next onset of menstruation and more daily temperature data are accumulated, the proposed method produced a considerably improved prediction (Figures (Figures33 and and4).4). During earlier stages of the cycle, however, previous methods tended to give better predictions. One of the reasons for this result is that we used the best possible prediction for each subject for the calendar calculation method as a baseline for comparison. Another reason is that, not being similar to the method proposed by Bortot et al. 6, the proposed method benefits little from the information about the lengths of previous cycles: the predictive distribution for upcoming menstruation does not vary largely (i.e., not being adjusted) at the days of onset of menstruation. Therefore, in practice, a combination of previous methods and the proposed method (e.g., using previous methods in earlier stages of the cycle, and then switching to the proposed method at some point in time) might result in a further improvement of prediction.

As measuring BBT is a simple and inexpensive means of determining the current phase of the menstrual cycle, the proposed framework can provide predictions of the next onset of menstruation that can easily be implemented for the management of women's health. In order to facilitate practitioners to implement and use the proposed method, we provide an r script for non‐Gaussian filtering and prediction, along with a simulated menstrual cycle dataset as a supporting web material. Although the non‐Gaussian filtering technique may be computationally impractical for a state‐space model with a high‐dimensional state vector 10, our proposed model only involves a two‐dimensional state vector so the computational cost required for filtering single data is negligible for a contemporary computer. We suggest that the proposed method will not be applicable to women using oral contraceptives because oral contraceptives can artificially control BBT and the onset of the menstruation. As long as the cycle is not under‐controlled by oral contraceptives and BBT shows the biphasic pattern, however, we suppose that the proposed prediction tool can be useful even for women with irregular cycle length (but see later for a difficulty in parameter estimation using data that include irregular cycle length).

The conditional probability density distribution for menstrual phase obtained with the sequential Bayesian filter may also be used for predicting other events that are related to the menstrual cycle. For example, the conditional probability of the menstrual phase being below or above a ‘change point’ of expected temperature could be used as a model‐based probability of the subject being in the follicular phase (before ovulation) or the luteal phase (after ovulation). These probabilities may also be used to estimate the timing of ovulation, for which BBT‐based estimates are error prone 11. Combined with a fecundability model that provides the probability of conception within a menstrual cycle that is conditional on daily intercourse behavior, the model might also be used to predict the probability of conception 6, 7, 8. As the Bayesian filtering algorithm can yield predictive and smoothed distributions (Equations (13) and (15)), these predictions could be obtained in both a prospective manner and a retrospective manner. Another possible extension of the model may be the addition of a new observation model for the physical condition of subjects, which may help to estimate the menstrual phase more precisely or to provide a forecast of the physical condition of subject as their menstrual cycle progresses.

We note that the state‐space model proposed in the present study may not be sufficiently flexible to describe the full variation in menstrual cycle length. In the present study, we fitted our model only to data from 20 subjects, and these data did not include extremely short or extremely long cycles because a preliminary investigation suggested that inclusion of such data could result in unreasonable parameter estimates (results not shown). Specifically, the estimates of α and β may become extremely small, leading to an almost flat probability density distribution for menstrual phase advancement. Under these parameter values, the predictive distribution for the onset of menstruation also becomes flat, preventing a useful prediction from being obtained. These results suggest that with a single set of parameter values, the proposed state‐space model does not capture the whole observed variation in cycle length. It is known that the statistical distribution of menstrual cycle length is characterized by a mixture distribution that comprises standard and nonstandard cycles, where, for the latter, a skewed distribution may well represent the observed pattern [e.g., 2, 3, 7, 8]. Our proposed model, however, does not provide such a mixture‐like marginal distribution for the day of onset of menstruation. It is also known that the mean and variance of the marginal distribution of cycle length can vary depending on subjects' age 3, 5, 6, 7. Therefore, we assume that modeling variations in the system model parameters (i.e., α and β), by including covariates such as age or within‐subject and among‐subject random effects, or both, may be a promising extension of the proposed framework. To this end, the model should be extended to treat longitudinal data (Y it,Z it), where time series is available for unit (i.e., subject or cycle) i. Then, differences in parameters among units can be modeled, for example, as logαi=γ0+jγjxij+ui, where x ij and u i are covariate (e.g., age) and random effect, respectively. To explain the skewed marginal distribution of menstrual cycle length, Sharpe and Nordheim 12 considered a rate process in which the development rate fluctuates randomly. Inclusion of random effects for the system model parameters would allow us to accommodate this idea. Although the inclusion of random effects will increase the flexibility of the framework, it would also considerably complicate the likelihood calculation and make parameter estimation more challenging.

The state‐space model proposed here provides a conceptual description of the menstrual cycle. We explain this perspective by using the analogy of a clock that makes one complete revolution in each menstrual cycle. The system model expresses the hand of the clock (i.e., the latent phase variable) moving forward with an almost steadily, yet slightly variable, pace. Observation models describe observable events that are imprinted on the clock's dial; for example, menstruation is scheduled to occur when the hand of this clock arrives at a specific point on the dial. The fluctuation in BBT, as well as other possibly related phenomena such as ovulation, would also be marked on the dial. Even though the position of the hand of the clock is unobservable, we can estimate it as a conditional distribution of the latent phase variable by using the Bayesian filtering technique, which enables us to make a sequential, model‐based prediction of menstruation. Although this is a rather phenomenological view of the menstrual cycle, it is useful for developing a rigorous and extendable modeling framework for predicting and studying phenomena that are associated with the menstrual cycle.

Supporting information

Supporting info item

Acknowledgements

We are grateful to K. Shimizu, T. Matsui, A. Tamamori, M. L. Taper, and J. M. Ponciano who provided valuable comments on this research. This manuscript benefited from comments and suggestions provided by two anonymous reviewers and an associate editor. These research results were achieved by ‘Research and Development on Fundamental and Utilization Technologies for Social Big Data’, the Commissioned Research of National Institute of Information and Communications Technology (NICT), Japan (no. 178A02). This research was initiated by an application to the ISM Research Collaboration Start‐up program (no. RCSU2013‐08) from M. Kitazawa.

Appendix A. State estimation, likelihood calculation, and prediction using the non‐Gaussian filter

A.1. 

The statistical inferences described in Section 3 of the main text were based on an estimation of the conditional distribution of the latent menstrual phase. Although the conditional distribution can be obtained by using the recursive equations described in Section 3.2, the integrals in those formulae are generally not analytically tractable for non‐linear, non‐Gaussian, state‐space models. Therefore, to obtain approximate conditional distributions, filtering algorithms for general state‐space models are required. We use the non‐Gaussian filter proposed by Kitagawa 10, in which the continuous state space is discretized with equally spaced fixed points at which the conditional probability density is evaluated. With these discretized probability density functions, integrals are evaluated by using numerical integration. This filtering algorithm also yields an approximate log‐likelihood, making the maximum likelihood estimation of the model parameters possible. Detailed descriptions for this numerical procedure are provided, for example, in Kitagawa 10, 13 and de Valpine and Hastings 14.

To apply non‐Gaussian filtering to the state‐space model described in Section 3.1, it is computationally more convenient to consider the state space of a circular phase variable rather than a linear phase variable described in the main text, because the latter requires discretization of unbounded real space. Hence, in the following, we first reformulate the statistical modeling framework in terms of a circular latent phase variable. The reformulated model conceptually agrees with the original formulation, given the probability that the menstrual phase advancement per day being more than 1 (i.e., Pr(ϵ t>1)) is negligible. In practice, this condition is naturally attained in the estimation of the model parameters as long as unusually short menstrual cycles (e.g., 3 to 4 days) are not included in the dataset. Later, we provide a numerical procedure to obtain the conditional distributions, log‐likelihood, and the predictive distribution of the next day of onset of menstruation based on the non‐Gaussian filtering technique (note, all definitions and notations are the same as those used in Section 3 of the main text).

A1. Model with a circular state variable

Instead of using the real latent state variable for the menstrual phase, θtR, we reformulate the model with a circular latent variable of period 1, denoted by ω t[set membership][0,1). We consider that the correspondence between these two variables is characterized by the following projection: ω t=f(θ t)=θ tmod1. The system model for this circular variable is then expressed as

ωt=ωt1+ϵtmod1
(A.1)

ϵtGamma(α,β).
(A.2)

Under this assumption, given ω t−1, the conditional distribution of ω t is a wrapped gamma distribution with a probability density function:

p(ωt|ωt1)=wGamma(α,β)=βαΓ(α)expβδ(ωt,ωt1)Φexp(β),1α,δ(ωt,ωt1),
(A.3)

where δ(ωt,ωt1)=(ωtωt1)I(ωt>ωt1)(1+ωtωt1)I(ωtωt1) is the advancement of the menstrual phase on the circular scale and Φ(z,s,a)=k=0zk(a+k)s is Lerch's transcendental function 15.

The observation model for BBT, which is conditional on ω t, is expressed as

yt=a+m=1M(bmcos2mπωt+cmsin2mπωt)+et
(A.4)

etNormal(0,σ2),
(A.5)

leading to a normal conditional distribution of y t with a probability density function:

p(yt|ωt)=Normalμ(ωt),σ2=12πσ2expytμ(ωt)22σ2,
(A.6)

where μ(ωt)=a+m=1M(bmcos2mπωt+cmsin2mπωt).

Assuming that Pr(ϵ t>1) is negligible, we consider menstruation to have started when ω t ‘steps over’ 1. This is represented as

zt=0whenωt>ωt1
(A.7)

=1whenωtωt1.
(A.8)

We can write this deterministic allocation in a probabilistic manner that is conditional on (ω t,ω t−1), where z t follows a Bernoulli distribution:

p(zt|ωt,ωt1)=(1zt)I(ωt>ωt1)+ztI(ωtωt1).
(A.9)

The log‐likelihood of this model can then be expressed as

l(ξ;YT,ZT)=logp(y1,z1|ξ)+t=2Tlogp(yt,zt|Yt1,Zt1,ξ),
(A.10)

where

logp(y1,z1|ξ)=log0101p(y1|ω1)p(z1|ω1,ω0)p(ω1,ω0)dω1dω0,
(A.11)

and for t=2,…,T,

logp(yt,zt|Yt1,Zt1,ξ)=log0101p(yt|ωt)p(zt|ωt,ωt1)p(ωt,ωt1|Yt1,Zt1)dωtdωt1.
(A.12)

For t=1,…,T, the prediction, filtering, and smoothing equations for the aforementioned model are expressed as follows:

p(ωt,ωt1|Yt1,Zt1)=p(ωt|ωt1)p(ωt1|Yt1,Zt1)=p(ωt|ωt1)01p(ωt1,ωt2|Yt1,Zt1)dωt2
(A.13)

p(ωt,ωt1|Yt,Zt)=p(yt,zt|ωt,ωt1)p(ωt,ωt1|Yt1,Zt1)0101p(yt,zt|ωt,ωt1)p(ωt,ωt1|Yt1,Zt1)dωtdωt1=p(yt|ωt)p(zt|ωt,ωt1)p(ωt,ωt1|Yt1,Zt1)0101p(yt|ωt)p(zt|ωt,ωt1)p(ωt,ωt1|Yt1,Zt1)dωtdωt1
(A.14)

p(ωt,ωt1|YT,ZT)=p(ωt,ωt1|Yt,Zt)01p(ωt+1,ωt|YT,ZT)p(ωt+1|ωt)p(ωt+1,ωt|Yt,Zt)dωt+1=p(ωt,ωt1|Yt,Zt)01p(ωt+1,ωt|YT,ZT)dωt+1p(ωt|Yt,Zt)=p(ωt,ωt1|Yt,Zt)01p(ωt+1,ωt|YT,ZT)dωt+101p(ωt,ωt1|Yt,Zt)dωt1,
(A.15)

where for t=1, we set p(ω 1,ω 0|Y 0,Z 0) as p(ω 1,ω 0), which is the specified initial distribution for the phase.

We consider the conditional probability that the next menstruation has started by day k+t, given the phase state ω t, to be F(k|ωt)=PrΔk(t)1ωt. Then, F(k|ω t) represents the conditional distribution function for the onset of menstruation, and under the assumption of the state‐space model described earlier, this is given as

F(k|ωt)=1ωtg(x;kα,β)dx=1G(1ωt;kα,β).
(A.16)

The conditional probability function for the day of onset of menstruation, f(k|ω t), is then given as

f(k|ωt)=F(k|ωt)F(k1|ωt)=1G(1ωt;kα,β)1G1ωt;(k1)α,β=G1ωt;(k1)α,βG(1ωt;kα,β),
(A.17)

where we set F(0|ω t)=0.

The marginal distribution for the day of onset of menstruation, h(k|Y t,Z t), is then obtained with the marginal filtering distribution for the phase state, p(ω t|Y t,Z t), which is given as

h(k|Yt,Zt)=01f(k|ωt)p(ωt|Yt,Zt)dωt,
(A.18)

which is used to provide a point prediction for the day of onset of menstruation by finding the value of k that gives the highest probability, maxh(k|Yt,Zt).

A2. Numerical procedure for non‐Gaussian filtering

We let ω(i),i=1,…,N be equally spaced grid points with interval [0,1). We denote the approximate probability density of the predictive, filtering, and smoothed distributions evaluated at these grid points on the state space by p˜p(i,j,t),p˜f(i,j,t), and p˜s(i,j,t), which are defined respectively as

p˜p(i,j,t)=pωt=ω(i),ωt1=ω(j)|Yt1,Zt1
(A.19)

p˜f(i,j,t)=pωt=ω(i),ωt1=ω(j)|Yt,Zt
(A.20)

p˜s(i,j,t)=pωt=ω(i),ωt1=ω(j)|YT,ZT,
(A.21)

where for i,j=1,…,N and t=1,…,T. As the initial distribution, we also specify the probability density p˜p(i,j,0). We define the probability density function of the system model (Equation (A.3)) evaluated at each grid point as p˜m(i,j)=pωt=ω(i)|ωt1=ω(j). Then, the prediction, filtering, and smoothing equations, respectively, are expressed as

p˜p(i,j,t)=p˜m(i,j)p˜f(j,t1)
(A.22)

p˜f(i,j,t)=pyt|ωt=ω(i)pzt|ωt=ω(i),ωt1=ω(j)p˜p(i,j,t)1N2i=1Nj=1Npyt|ωt=ω(i)pzt|ωt=ω(i),ωt1=ω(j)p˜p(i,j,t)
(A.23)

p˜s(i,j,t)=p˜f(i,j,t)p˜s(i,t+1)/p˜f(i,t),
(A.24)

where p˜f(i,t) and p˜s(i,t+1) are the marginalized filtering and smoothed distributions, respectively, which are obtained as p˜f(i,t)=k=1Np˜f(i,k,t)/N and p˜s(i,t+1)=k=1Np˜s(k,i,t+1)/N, respectively. Note in practice that the predictive and smoothed distributions should be normalized so that the value of the integral over the whole interval becomes 1 13.

The denominator of Equation (A.23) provides the approximate likelihood for an observation at time t. Therefore, the log‐likelihood of the model is approximated as

l(ξ;YT,ZT)=t=1Tlogi=1Nj=1Npyt|ωt=ω(i)pzt|ωt=ω(i),ωt1=ω(j)p˜p(i,j,t1)2TlogN.
(A.25)

Finally, the marginal distribution for the day of initiation of menstruation is approximated as

h(k|Yt,Zt)=1Ni=1Nfk|ωt=ω(i)p˜f(i,t).
(A.26)

Appendix B. Details for comparing accuracy of the prediction for the day of onset of menstruation

B.1. 

We compared the root mean square error (RMSE) and the mean absolute error (MAE) of prediction of the next day of onset of menstruation among the proposed method, calendar calculation method, and the method proposed in Bortot et al. 6.

As described in the main text, the proposed method can adaptively adjust the prediction using the within‐cycle information (i.e., daily BBT records). Thus, we estimated RMSE and MAE of the proposed method for several points in time within the cycle. Specifically, we calculated RMSE and MAE of the prediction obtained at the day of onset of the previous menstruation, as well as those obtained at 21, 14, 7, 6, 5, 4, 3, 2, and 1 day(s) before the upcoming day of onset of menstruation.

In the calendar calculation method, RMSE and MAE depend on the fixed number of days used for prediction. Thus, for each subject, we calculated RMSE and MAE over a range of parameter values and used the lowest value as a baseline for comparison (Tables 3 and 4).

The method proposed by Bortot et al. 6 uses state‐space modeling framework to predict the length of the cycle. The predictive distribution of menstrual cycle length is obtained conditional on past time series of cycle length. In contrast to the proposed method, the prediction does not vary within the cycle.

In order to implement the prediction method of Bortot et al. 6, we first fitted state‐space models to data for parameter estimation (Table 1). Samples from the joint posterior distribution of parameters and state variables were obtained by using the Markov chain Monte Carlo (MCMC) method run on jags software 16. Then, the ensemble Kalman filter was applied to obtain one‐step‐ahead predictive distributions of cycle length. The MCMC outputs were used for parameters and initial distributions of state variables that are required for the ensemble Kalman filter to work.

We tested several model variants that were simplified versions of the original state‐space model described in Bortot et al. 6. Among model variants in which convergent MCMC samples were obtained, we identified a model that provides, on average, the most accurate predictions and used it as a baseline for the comparison (Figures (Figures33 and and44).

Notes

Fukaya K., Kawamori A., Osada Y., Kitazawa M., and Ishiguro M. (2017) The forecasting of menstruation based on a state‐space modeling of basal body temperature time series. Statist. Med., 36: 3361–3379. doi: 10.1002/sim.7345.

References

1. Barron ML, Fehring RJ. Basal body temperature assessment: is it useful to couples seeking pregnancy? American Journal of Maternal Child Nursing 2005; 30(5):290–296. [PubMed]
2. Harlow SD, Zeger SL. An application of longitudinal methods to the analysis of menstrual diary data. Journal of Clinical Epidemiology 1991; 44(10):1015–1025. [PubMed]
3. Guo Y, Manatunga AK, Chen S, Marcus M. Modeling menstrual cycle length using a mixture distribution. Biostatistics 2006; 7(1):100–114. [PubMed]
4. Lin X, Raz J, Harlow SD. Linear mixed models with heterogeneous within‐cluster variances. Biometrics 1997; 53(3):910–923. [PubMed]
5. Huang X, Elliott MR, Harlow SD. Modelling menstrual cycle length and variability at the approach of menopause by using hierarchical change point models. Journal of the Royal Statistical Society: Series C (Applied Statistics) 2014; 63(3):445–466. [PMC free article] [PubMed]
6. Bortot P, Masarotto G, Scarpa B. Sequential predictions of menstrual cycle lengths. Biostatistics 2010; 11(4):741–755. [PubMed]
7. McLain AC, Lum KJ, Sundaram R. A joint mixed effects dispersion model for menstrual cycle length and time‐to‐pregnancy. Biometrics 2012; 68(2):648–656. [PubMed]
8. Lum KJ, Sundaram R, Louis GMB, Louis TA. A Bayesian joint model of menstrual cycle length and fecundity. Biometrics 2016; 72(1):193–203. [PubMed]
9. Särkkä S. Bayesian Filtering and Smoothing. Cambridge University Press: Cambridge, 2013.
10. Kitagawa G. Non‐Gaussian state‐space modeling of nonstationary time series. Journal of the American Statistical Association 1987; 82(400):1032–1041.
11. Dunson DB, Weinberg CR. Modeling human fertility in the presence of measurement error. Biometrics 2000; 56(1):288–292. [PubMed]
12. Sharpe PJH, Nordheim AW. Distribution model of human ovulatory cycles. Journal of Theoretical Biology 1980; 83(4):663–673. [PubMed]
13. Kitagawa G. Introduction to Time Series Modeling. Chapman & Hall/CRC: Boca Raton, 2010.
14. de Valpine P, Hastings A. Fitting population models incorporating process noise and observation error. Ecological Monographs 2002; 72(1):57–76.
15. Coelho CA. The wrapped Gamma distribution and wrapped sums and linear combinations of independent Gamma and Laplace distributions. Journal of Statistical Theory and Practice 2007; 1(1):1–29.
16. Plummer M. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. Proceedings of the 3rd international workshop on distributed statistical computing (dsc 2003), Vol. 124. Technische Universit at Wien, Austria, 2003, 125.

Articles from Wiley-Blackwell Online Open are provided here courtesy of Wiley-Blackwell, John Wiley & Sons