Biostatistics Advance Access originally published online on April 14, 2005
Biostatistics 2005 6(3):434-449; doi:10.1093/biostatistics/kxi020
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
A transformation approach for incorporating monotone or unimodal constraints
Jiann-Ping Hsu School of Public Health, Georgia Southern University, P.O. Box 8076, Statesboro, GA 30460-8076, USA lgunn{at}georgiasouthern.edu
Biostatistics Branch, National Institute of Environmental Health Sciences, MD A3-03, P.O. Box 12233, Research Triangle Park, NC 27709, USA
* To whom correspondence should be addressed.
 |
SUMMARY |
|---|
 TOP SUMMARY 1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
Samples of curves are collected in many applications, including studies of reproductive hormone levels in the menstrual cycle. Many approaches have been proposed for correlated functional data of this type, including smoothing spline methods and other flexible parametric modeling strategies. In many cases, the underlying biological processes involved restrict the curve to follow a particular shape. For example, progesterone levels in healthy women increase during the menstrual cycle to a peak achieved at random location with decreases thereafter. Reproductive epidemiologists are interested in studying the distribution of the peak and the trajectory for women in different groups. Motivated by this application, we propose a simple approach for restricting each woman's mean trajectory to follow an umbrella shape. An unconstrained hierarchical Bayesian model is used to characterize the data, and draws from the posterior distribution obtained using a Gibbs sampler are then mapped to the constrained space. Inferences are based on the resulting quasi-posterior distribution for the peak and individual woman trajectories. The methods are applied to a study comparing progesterone trajectories for conception and nonconception cycles.
Keywords: Changepoint; Gibbs sampler; Hormone measurements; Menstrual cycle; Nested data; Order-restricted inference; Progesterone; Shape constraint
|
1. INTRODUCTION |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
Reproductive epidemiologists and endocrinologists are interested in studying hormone patterns in the menstrual cycle. Differences in hormone levels among women during the various phases of the cycle may be indicative of reproductive function (Baird et al., 1997
, 1999). For certain hormones, such as progesterone, studies show that the trajectories increase monotonically to a peak at an unknown location with decreases thereafter (Yen and Jaffe, 1986). Characteristics of these curves, including peak location, may provide insight into the probabilities of conception and early pregnancy loss (Baird et al., 1997, 1999).
In such biomedical studies, smoothing methods are commonly used to estimate mean trajectories, while limiting parametric assumptions. Zhang et al. (1998, 2000) used semiparametric models for correlated data, and estimated smoothing parameters using restricted maximum likelihood estimation. Brumback and Rice (1998) used smoothing splines with mixed effects models to analyze progesterone data in the menstrual cycle. Although such approaches have many advantages, they did not incorporate shape constraints on the mean curves. Such constraints provide a means for incorporating biological information into the analysis without introducing restrictive parametric assumptions about the unknown trajectories. Restricting the shapes that are possible to a class consistent with prior information can improve statistical efficiency while producing a more realistic estimate.
Several approaches have been proposed for estimating curves subject to shape constraints. Gelfand and Kuo (1991) and Ramgopal et al. (1993) used restricted nonparametric priors to accommodate shape restrictions in potency curves. Gelfand et al. (1992) proposed a general Gibbs sampling strategy for incorporating parameter constraints by discarding draws inconsistent with the constraint. Lavine and Mockus (1995) proposed a general method for nonparametric Bayesian isotonic regression. Dunson and Columbo (2003) proposed a hierarchical model and smoothing prior that allowed for shape restrictions on mean curves underlying categorical observations. Adaptation of these approaches to the problem of incorporating constraints on higher level parameters in a hierarchical model (e.g. individual-specific means) is not straightforward.
In addition to these Bayesian approaches, frequentist methods are available for monotone curve estimation. Mammen (1991) used smoothing and isotonization steps to efficiently estimate a monotone regression function. Mammen and Thomas-Agnan (1999) proposed a method for constraining smoothing splines by first calculating the unrestricted smoothing spline, and then projecting onto the constrained space using a Sobolev-type norm. Such isotonization tends to yield a smaller mean square error in many cases (Robertson et al., 1988; Mammen, 1991; Hwang and Peddada, 1994).
Applying the isotonization idea to the Bayesian paradigm, Dunson and Neelon (2003) proposed an approach for order-restricted inference in generalized linear models. Instead of mapping a single unrestricted estimate to the constrained space, they proposed applying an isotonic regression transformation to draws from the unconstrained posterior density obtained using a Markov chain Monte Carlo (MCMC) algorithm. Their focus was on simple ordering in regression parameters, and not on ordering with unknown changepoints or constraints on higher level parameters in a hierarchical model.
Motivated by applications to hormone studies, we generalize the method of Dunson and Neelon (2003) to allow for umbrella ordering and constraining of higher level parameters in a hierarchical model. We use the term ‘umbrella’ to refer to a monotone or unimodal ordering, since it is widely used in the literature on order-restricted inference. To generalize isotonic regression transformations to the case where the peak location is unknown, we choose the peak that minimizes the weighted least squares distance between the unconstrained and umbrella-ordered parameters. We then apply this transformation to individual-specific means generated from a hierarchical model.
Section 2 considers umbrella-ordered cases when data are independent and normal. Section 3 generalizes the approach to constrain individual-specific trajectories characterized by a hierarchical model. Section 4 applies the method to progesterone data considered previously by Brumback and Rice (1998) and provides a simulation example, while Section 5 discusses the results.
|
2. UMBRELLA-ORDER RESTRICTIONS |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
We focus initially on the inference for a vector of normal means known a priori to follow an umbrella-order restriction. Let 
where 
is the mean vector and 
= 
2Ip x p is the covariance matrix, for
. Ignoring the ordering initially, a conjugate prior for µ can be specified as 
where 
The posterior of µj conditional on 2 and y is then 
where 
and 
We focus inference on an order-constrained parameter µ*, which arises from a mapping of the unconstrained mean vector µ from 
where 
is a subset defined by the following set of inequalities on the elements of µ such that 
We define 
as the set of umbrella-ordered means with peak 
:
 | (2.1) |
When the peak occurs at 1 or p, a simple nonincreasing or nondecreasing order, respectively, results. The support of µ* is therefore the space of p x 1 vectors following an umbrella ordering with peak location
Following Dunson and Neelon (2003) who considered simple ordering, we choose a mapping from 
p 
that minimizes the distance between µ* and µ. In particular, based on Robertson et al. (1988) and Hwang and Peddada (1994), we use an isotonic regression transformation that sets the elements of µ* equal to weighted averages of the elements of µ. Assuming a peak at , the transformation follows the form shown in Robertson et al. (1988) and Dunson and Neelon (2003):
 | (2.2) |
where µ|, y is the unconstrained posterior covariance and 
and 
denote subsets of 
such that the ordering µj' 
µj is known for all 
and the ordering µj' 
µj is known for all 
The [s:t] subscript represents the submatrices and subvectors corresponding to elements s, s + 1, ..., t. Note that in this simple case, 
To allow for a peak at an unknown location, , we choose µ* by minimizing the Mahalanobis distance measure across different choices of peak:
 | (2.3) |
Using this expression in combination with expression (2.2) produces the value of µ* 
, which is as close as possible to the unrestricted value 
The weighted least squares distance measure minimizes the distance between µ* and µ in Euclidean space, requiring the distance to be relatively small for elements of µ about which much is known (i.e. the posterior variance is small). In this manner, the mapping tends to result in values of µ* with high density in the posterior for µ. We follow Dunson and Neelon (2003) in taking the view that µ* is an order-restricted functional of the parameter µ, and hence, we can consider the draws of µ* as draws from a Bayesian posterior.
An alternative approach would be to directly place a prior distribution and define a likelihood based on µ*. Although this more conventional Bayesian strategy is conceptually appealing, our transformation approach has important practical advantages. In particular, it is typically the case that investigators are interested in both the order-restricted µ* and the unrestricted µ, since assessing changes in the estimates due to the incorporation of the constraint is a critical component of the analysis (particularly, in applications in which one may be somewhat skeptical about the constraint). In addition to computational advantages in more complex problems (as described in Sections 3 and 4) and a tendency to produce estimates that are close to unrestricted ridge-type estimators, the transformation strategy has the advantage that it is straightforward to assess the impact of incorporating the constraint. One can simply examine how inferences change using µ* instead of µ.
The approach is also appealing from a frequentist perspective due to the close relationship with traditional isotonic regression transformations. In addition, one can show that the mode of the posterior for µ* converges to the constrained estimator in the limit as n 
. Using this equivalence, it is also straightforward to show that the mode of the posterior for µ* is a consistent estimator as long as the true means belong to . This proof follows along the lines described by Turner and Wollan (1997) for their classical estimator, and is one motivation for our use of the Mahalanobis distance measure in choosing the peak. A final comment is that the posterior tends to be centered close to the constrained estimator when the information in the prior is small relative to the sample size even when the sample size is small to moderate. In complex settings in which it is very difficult to estimate the constrained estimator, the transformation strategy can be viewed as a reasonable alternative.
For additional theoretical properties, we focus on the simple p = 3 case, which is not unusual to find in dose-response and bioassay applications. Theorems for general p are very difficult to establish (as in classical order-restricted inference problems).
THEOREM 1 For p(see Appendix for the proof).
This theorem implies that the transformation from 
maintains the peak location. Since we wish to limit systematic biases in mapping from this is an appealing property.
We compare the proposed transformation approach from Section 2.1 with the standard Bayesian approach of placing a truncated conjugate prior on the restricted space. A typical prior used to constrain µ follows the truncated normal form:
 | (2.4) |
where µ0 and 0 are defined in Section 2.1. Conditional on 
and the peak location occurring at , the posterior distribution for µ arising from prior (2.4) follows the truncated normal distribution:
 | (2.5) |
where f(·) and F(·) are the normal density and distribution functions, respectively, and 
and 
are defined in Section 2.1, with 
when the prior variance is large. Sampling from (2.5) is equivalent to drawing from the unrestricted posterior 
and discarding values that are less than max{µ–}. This conditional posterior gives zero prior probability to µ max{µ–} and has expectation
 | (2.6) |
As max{µ–} increases and the conditional posterior variance decreases, the expectation in (2.6) can be substantially larger than and hence, when the prior variance is large, substantially larger than 
the mean of the data.
It is of interest to compare expression (2.5) with the conditional posterior distribution for 
obtained under the proposed transformation approach:
 | (2.7) |
which is equivalent to sampling from the unrestricted 
posterior and setting values less than max{µ–*} equal to the boundary value max{µ–*}. It follows that the conditional estimator using the transformation approach is
 | (2.8) |
which is closer to the unrestricted posterior mean, 
than the conditional posterior mean shown in expression (2.6) for the truncated normal prior. This result suggests that when the prior precision, 
is small relative to the information in the data, the proposed transformation approach yields estimates closer to the data and less subject to systematic biases compared with the more conventional Bayesian approach. This difference in the two estimates can be substantial, particularly when the true value is close to the boundary of , when the sample size is small to moderate, and when strong constraints are assigned. Note that we have focused on the case where the peak is known. When the peak location is unknown, the standard strategy would be to place a prior distribution (e.g. uniform) on the peak location and estimate this changepoint simultaneously with the mean values. It is well known that the introduction of an unknown changepoint can greatly contribute to the computational burden. This is not a major hurdle in the simple ordered normal means case but is in cases where the changepoint can vary for each subject, as in the menstrual cycle application.
It is useful to consider a simple illustrative example. Suppose p = 3 and that the group-specific empirical means are 
with all three groups having equal sample sizes. We fix the empirical variance of y at 100. It is known a priori that there exists an umbrella ordering in µ with the peak location unknown. Figure 1 plots the estimated posterior means for µ1, µ2, and µ3 for a range of values of c and n under the approach which assigns a uniform improper prior for µ with support on and plots the estimators under our proposed transformation approach. As expected, the first and third panels are identical due to symmetry.
|
It is apparent from Figure 1 that the estimators under the transformation approach are closer to the empirical means in each of the cases, with the largest differences occurring for small sample sizes and negative values of c. Since having negative values for c produces a bathtub-shaped ordering, the umbrella-ordered constraint is clearly violated. Figure 1 suggests that, even in the case where the prior is chosen to be noninformative subject to the order restriction and the data do not violate the order constraint, the posterior densities produced by the standard approach can be centered on values far from the empirical means of the data. The transformation approach is far less sensitive to this problem. Although we have focused on small samples in this illustrative example, this problem is also important in large sample settings in which we are interested in estimating individual-specific means but do not have a large amount of information on each individual.
|
3. ORDER RESTRICTIONS IN HIERARCHICAL MODELS |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
Until this point, we have focused on the simple case where the data consist of independent observations from a normal distribution. However, it is straightforward to apply the methods to any hierarchical model in which umbrella restrictions are needed. We focus on the case where there are repeated multivariate observations on each of n subjects, so that outcome data on subject i consist of a vector of observations Yij = (yij1, ... yijp)' at times j = 1, ..., ni. For example, in the menstrual cycle application, i indexes woman and j indexes menstrual cycle within a woman. Let µij = (µij1, ..., µijp)' denote the linear predictor in an arbitrary hierarchical generalized linear mixed model. This linear predictor is characterized as a function of fixed effects regression parameters,
, individual and potentially time-specific random effects, ßi, variance component parameters, 
, and residual error variances, 
. A prior distribution for µij is induced by choosing a random effects density and prior distributions for the population parameters. The usual Bayesian approach of placing constraints on the population parameters by explicitly choosing priors with restricted support is straightforward, with some risk of bias (as discussed in Section 2). However, restricting the population parameters is not the same as restricting the linear predictor µij , which is our focus.
It is difficult to induce restrictions on µij through restricting the population parameters and random effects, since constraints on the higher level parameters imply complex restrictions on the lower level parameters. One can potentially use an MCMC algorithm in which, for each parameter one at a time (including the random effects), the appropriate constraint is calculated conditionally on current values of the other parameters, and the parameter is then updated subject to this constraint. This approach has several important problems. First, the parameter constraint is often difficult to calculate being dependent on a function of the other parameters and the predictors. This calculation is made more challenging if one wants to restrict the parameters so that µij for values of the predictors not observed in the study sample. Second, since the constraint on a particular parameter is highly dependent on the current values of the other parameters, there will tend to be high levels of autocorrelation in the Markov chain and hence poor computational efficiency. In addition, we do not want to restrict the peak in µij to occur at the same location for all i, j. The standard Bayesian approach requires the explicit specification of a hierarchical prior for the peak location. In conducting posterior computation, it is then necessary to update the parameters characterizing the distribution of the peak along with the individual-specific changepoint parameters. Since efficient computation is difficult even for models with a single changepoint, computation for umbrella-ordered curves with varying peaks is extremely challenging using the standard methods.
The transformation approach proposed in Section 2 effectively solves these problems. We begin by collecting draws from the unconstrained posterior distribution of 
= (, ß, , ) using a standard Gibbs sampling algorithm without incorporating the constraint. This is straightforward, and the parameters can be updated in blocks to improve computational efficiency. The linear predictor µij is then computed as a function of for each draw collected from the unconstrained posterior. Once µij is generated, the transformation approach is used to map µij µij* .
|
4. APPLICATION TO MENSTRUAL CYCLE DATA |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
We illustrate the methods through application to data from a study of progesterone levels in the menstrual cycle, analyzed previously by Brumback and Rice (1998). The data consist of progesterone metabolite measurements in daily urine samples of women with healthy reproductive function. These hormone measurements are recorded on the log scale as a variance stabilizing mechanism. Days are nested within cycles which are nested within women who are in one of the two groups—nonconceptive or conceptive. Of the 51 women in the study, 29 of them belong to the nonconceptive group (69 cycles) and 22 belong to the conceptive group (22 cycles). Among the nonconceptive group, women have between one and five cycles each, with 2 as the median number of cycles. Data are aligned relative to the day of ovulation estimated using urinary luteinizing hormone (LH) surge. We focus on the 24 day window starting 8 days prior to ovulation and ending 15 days after ovulation, since interest focuses on comparing conceptive and nonconceptive cycles with respect to pre- and postovulatory progesterone patterns. Ovulation is referred to as day 0. Notice that we are not standardizing cycles to 28 days as is commonly done throughout the literature, because that results in a loss of information and does not seem defensible biologically.
Figure 2 shows the hormone measurements for two nonconceptive women and for two conceptive women. As expected, progesterone levels tend to increase to a peak and then decrease in nonconceptive cycles, while in conceptive cycles increases continue throughout the window. Note that both of these patterns are consistent with an umbrella ordering, with the peak occurring at the end of the window for conceptive cycles and earlier for nonconceptive cycles.
|
Let yijk represent the hormone measurement for day k = {1, ..., 24} of cycle j = {1, ..., ni} from woman i = {1, ..., 51}. We first choose a simple autoregressive hierarchical model to characterize the data without considering the constraint. We focus on a discrete time model, since hormone measurements are collected daily. Given the possible diurnal variability and standardization of collection protocols to focus on first morning urine samples, it is not appropriate to extrapolate to a finer scale based on the daily observations.
 | (4.1) |
where cij is an indicator of the conception group for the jth cycle of the ith woman, 
are day-specific population means [in either the conception (C) or the nonconception (N) group], 
are women-specific deviations from the population mean, 
measures residual autocorrelation, and 
are the residual errors. To complete the prior specification, 
is the woman-specific intercept, 
measures heterogeneity in the intercept, and 
measures autocorrelation in the ßik. The inverse gamma prior parameters for 2, 
and 
are chosen to be moderately diffuse with a mean of 5 and variance of 100, and a sensitivity analysis for different prior parameter choices was conducted. The prior mean of , 0, is set equal to the empirical mean from a previous study conducted by Baird et al. (1999). The prior variance for is set at 10 to allow for differences in the studies.
Also notice the choice of prior density for ßik. While it is common to place a normal prior centered at 0 with a huge variance on the ßik values, this approach gives high probability to values that are known to be implausible, resulting in slow mixing of the Gibbs sampler. For this reason, we choose the prior specification for ßik as mentioned above. Furthermore, ßik is chosen over ßi or ßij because ßik allows each woman to have a distinct pattern in her day-specific means.
Our focus lies in the order-restricted functional, µi* = gi(µi), of the unconstrained women-specific means, where: 
The ordering in 
parallels that of published data and results (Baird et al., 1997, 1999). The restriction on guarantees that the progesterone hormone measurement in a woman's menstrual cycle increases monotonically to an unknown changepoint, that varies from woman to woman, with decreases thereafter.
Since the full conditional distributions of each of the unconstrained parameters follow a simple conjugate form, posterior computation can proceed via a simple Gibbs sampling algorithm which proceeds as follows:
Step 1:
Missing observations are sampled from their full conditional posterior distribution under the missing at random assumption.
Step 2:
Sample
Step 3:
Calculate women-specific means, µik, as a deterministic function of the day-specific population means and the women-specific random effects.
Step 4:
The unconstrained µik values are mapped to the umbrella-ordered space
The algorithm was run for 110 000 iterations. The first 10 000 iterations were discarded as a burn-in, and every 100th iteration of the following 100 000 samples were saved to thin the chain. Traceplots indicate convergence of the algorithm and good mixing.
We show unconstrained and constrained mean estimates and 95% credible intervals for two conceptive and two nonconceptive women in Figure 2. Notice in Figure 2 and Table 1 that the unconstrained estimates have bumps which likely reflect noise in the assay and not true changes in the woman-specific means. Table 1 shows day-specific means for women in each group. The constrained estimates do not vary systematically from the unconstrained estimates, implying minimal bias in the estimates. As expected, the interval estimates under the transformation approach are generally narrower than the unconstrained credible intervals, particularly prior to ovulation, reflecting a possible improvement in efficiency for the constrained approach. These differences in the unconstrained and constrained point and interval estimates are also seen in the simulation study in Section 4.5. Furthermore, a likely improvement in efficiency for the constrained approach is also illustrated in Figure 3, which shows the ratios of the unconstrained to constrained estimated posterior variances in the women-specific curves. In particular, we estimate a hormone curve for each woman, and there is uncertainty in estimating these women-specific curves. Figure 3 quantifies how this uncertainty decreases under the order-restricted procedure. As with the narrower interval widths displayed in Figure 2 corresponding to the constrained means, Figure 3 shows a likely improvement in efficiency, particularly prior to ovulation.
|
|
Epidemiologists are also interested in the peak location and in systematic differences in this peak between nonconceptive and conceptive cycles. The day at which progesterone levels peak in nonconceptive cycles is thought to be close to the day of implantation in conceptive cycles. The timing of the peak may be informative about the probability of conception, and women who conceive after progesterone levels fall are at a greater risk of early pregnancy loss and later adverse outcomes (Baird et al., 1999
). For this reason, trajectory characteristics, including peak location, are very interesting. Figure 4 provides a histogram of peak locations for nonconceptive and conceptive cycles. Since left-skewness in the posterior distribution caused the posterior mean estimator of the peak to occur prior to the end of the cycle in conceptive women, the posterior mode is shown in Figure 4. The most likely peak day occurs 5 days after ovulation for the nonconceptive women and 15 days after ovulation for the conceptive women. The peak for the conceptive women occurs at the end of the cycle window. Furthermore, there is an approximate posterior probability of 1.0 that the peak occurs later on average for conceptive cycles than for nonconceptive cycles.
|
Another interest of epidemiologists is comparing conceptive trajectories with nonconceptive trajectories during certain days of the cycle. Breaking the cycle days into three intervals: (1) [– 8, 1], (2) [2, 8], and (3) [9, 15] (with day 0 corresponding to ovulation day), we find that the estimated probabilities of higher mean progesterone levels for conception cycles are 0.46, 0.01, and 0.69 for intervals 1–3, respectively. Therefore, conceptive and nonconceptive cycles are similar except during the week following ovulation during which conceptive cycles have significantly lower progesterone levels. There is no evidence of differences in preovulatory progesterone levels between conceptive and nonconceptive cycles, suggesting that progesterone levels during the fertile interval have little effect on fecundability. Note that we have ended the first interval the day after the estimated day of ovulation, which was based on urinary LH surge in this study, to allow for possible measurement error in the LH estimate.
Aside from characteristics of mean trajectories, estimates for 2, 1, and are also computed. The posterior mean for the residual variance is 
with 95% credible interval (0.23, 0.26). The autocorrelation parameter, 1, has a posterior mean of 0.65, with 95% credible interval (0.60, 0.70). The heterogeneity in the intercept, and in the women-specific random effects, 
have posterior means of 0.95 and 0.10, with 95% credible intervals (0.56, 1.50) and (0.08, 0.12), respectively.
Finally, we conducted a sensitivity analysis to assess the robustness of these results to the prior specification. We varied the prior means for 2, and from 1 to 7 and the prior variances from 10 to 10 000, allowing arbitrarily smaller and larger values than those used in the primary analysis. Overall, posterior estimates for these variance component parameters were robust to the prior specification, with the estimate for 
being the most sensitive. The parameters of primary interest, the restricted woman-specific means, varied approximately 1% on average across the different prior specifications.
In order to check our computational algorithm and methods as well as explore frequentist operating characteristics, we ran the analysis using 25 simulated data sets having the same structure as the progesterone application data but with simulated progesterone values. The constrained means from the progesterone analysis were used as the true values for the simulated data. In implementing the analysis, we used a burn-in period of 10 000 iterations, with the chain run an additional 50 000 iterations, collecting samples at every 100th iteration.
Based on the results for the 25 simulated data sets, we calculated pointwise bias, defined as the average difference between the restricted estimates and the true values of the parameters, in the constrained day-specific mean estimates as well as the estimate for peak location, and the results are shown in Tables 2 and 3. There is no evidence of systematic bias in the estimates, and the differences between the averaged restricted means and the true values may reflect Monte Carlo error. Also notice the minimal bias in peak location across women and analyses in the two groups, with all biases less than 0.75 days. Table 2 also shows the average pointwise 95% credible intervals for the restricted day-specific means across the 25 analyses, and coverage. Notice the differences in slightly lower and higher coverage prior to ovulation and postovulation, respectively. It appears that coverage of 95% credible intervals is consistent with the nominal level, though a larger simulation study is needed to fully assess this.
|
|
|
5. DISCUSSION |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
This article introduces a useful and easy-to-implement method for constraining parameters to follow an umbrella ordering. In the progesterone application, each cycle's hormone curve can be viewed as a single random function, with repeated random functions then collected for each woman. This is an example of hierarchical functional data, which refer to correlated samples of functions, as opposed to standard hierarchical data, which refer to data having some natural dependency structure, such as cycles within women. Although we have applied this approach to hierarchical functional data, there are many other settings in which the methodology is potentially useful.
The transformation approach can be used for inferences on dose-response functions with possible downturns at higher dose levels, or hazard functions with bathtub or hill shapes. Applying the approach in other settings is trivial since one can save existing output from an unconstrained analysis (using WinBUGS, for example), and transform the output using a simple function (i.e. coded in S-Plus or SAS). In contrast, following more conventional Bayesian procedures of placing an explicit prior with constrained support requires new coding for each new case, which is highly challenging to implement for hierarchical models with umbrella constraints on higher level parameters or for dose-response functions mentioned above.
We use an empirical Bayes-type procedure based on transforming draws from an unconstrained posterior distribution for the parameters. Although this procedure may correspond to a fully Bayesian procedure for a particular choice of implicit prior distribution for the parameters, such a prior is extremely challenging to derive, except in very simple, special cases. Hence, the proposed approach should not be considered as a fully Bayesian procedure.
Since computational hurdles are one of the major factors preventing widespread use of order and shape restrictions in analyses of biomedical data, our method is all the more appealing for its promise of computational efficiency. In addition to computational advantages, our results suggest that we can limit bias induced by the order restriction by using a mapping which minimizes the distance between the restricted and unrestricted estimates. Placing order constraints on parameters is a useful way in which outside information can be incorporated into the analysis in a nonparametric manner without using a highly informative prior for the specific values of the parameters. Such constraints often result in smoothing of the estimates and improvements in efficiency.
|
APPENDIX |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
For p = 3, there are six possible inequalities in the means:
 |
In general, the values of the constrained parameters given a peak at = 1 are defined as follows:
 |
 |
 |
where hl are the elements of the unconstrained posterior covariance, µ|, y, for h, l = 1, ..., p. Next, µ*1 is computed using the above results for each of the six inequalities, respectively:
 |
 |
 |
 |
 |
 |
This procedure allows us to find the inequality for which the Mahalanobis distance (2.3) is minimized, hence providing an ordering in the restricted estimates, µ*. For the ordering in which the peak occurs at 1, this process shows the loss corresponding to the second inequality, µ1 µ2 µ3, is 0. No other loss pertaining to the remaining inequalities under the ordering µ1 µ2 µ3 in the unconstrained estimates has a value this small. Comparing the distances to find the one that satisfies expression (2.3) is not always as simple as with this ordering in the unconstrained estimates specified here. When the ordering in the unrestricted estimates is different from the one specified here, not all of the inequalities yielding the smallest value of (2.3) attains a loss of 0. Nonetheless, computing (2.3) for each of the orderings in the unconstrained estimates will show that the ordering in µ* is equivalent to the ordering in µ for p 3.
|
ACKNOWLEDGMENTS |
|---|
We thank Donna Baird, Chong Tu, and Dalene Stangl for useful comments and discussions. We also thank Babette Brumback for providing data for the example.
|
REFERENCES |
|---|
|
TOPSUMMARY1. INTRODUCTION2. UMBRELLA-ORDER RESTRICTIONS3. ORDER RESTRICTIONS IN...4. APPLICATION TO MENSTRUAL...5. DISCUSSIONAPPENDIXREFERENCES |
|---|
-
BAIRD, D. D., WEINBERG, C. R., ZHOU, H., KAMEL, F., MCCONNAUGHEY, D. R., KESNER, J. S. AND WILCOX, A. J. (1999). Preimplantation urinary hormone profiles and the probability of conception in healthy women. Fertility and Sterility 71, 40–49.[CrossRef][Web of Science][Medline]
BAIRD, D. D., WILCOX, A. J., WEINBERG, C. R., KAMEL, F., MCCONNAUGHEY, D. R., MUSEY, P. I. AND COLLINS, D. C. (1997). Preimplantation hormonal differences between the conception and non-conception menstrual cycles of 32 normal women. Human Reproduction 12, 2607–2613.
BRUMBACK, B. A. AND RICE, J. A. (1998). Smoothing spline models for the analysis of nested and crossed samples of curves. Journal of the American Statistical Association 93, 961–976.[CrossRef][Web of Science]
DUNSON, D. B. AND COLOMBO, B. (2003). Bayesian modeling of markers of day-specific fertility. Journal of the American Statistical Association 98, 28–37.[CrossRef]
DUNSON, D. B. AND NEELON, B. (2003). Bayesian inference on order-constrained parameters in generalized linear models. Biometrics 59, 286–295.[CrossRef][Web of Science][Medline]
GELFAND, A. E. AND KUO, L. (1991). Nonparametric Bayesian bioassay including ordered polytomous response. Biometrika 78, 657–666.
GELFAND, A. E., SMITH, A. F. M. AND LEE, T. M. (1992). Bayesian analysis of constrained parameter and truncated data problems using Gibbs sampling. Journal of the American Statistical Association 87, 523–532.[CrossRef]
HWANG, J. T. G. AND PEDDADA, S. D. (1994). Confidence interval estimation subject to order restrictions. Annals of Statistics 22, 67–93.
LAVINE, M. AND MOCKUS, A. (1995). A nonparametric Bayes method for isotonic regression. Journal of Statistical Planning and Inference 46, 235–248.[CrossRef]
MAMMEN, E. (1991). Estimating a smooth monotone regression function. The Annals of Statistics 19, 724–740.
MAMMEN, E. AND THOMAS-AGNAN, C. (1999). Smoothing splines and shape restrictions. Scandinavian Journal of Statistics 26, 239–252.
RAMGOPAL, P., LAUD, P. W. AND SMITH, A. F. M. (1993). Nonparametric Bayesian bioassay with prior constraints on the shape of the potency curve. Biometrika 80, 489–498.
ROBERTSON, T., WRIGHT, F. AND DYKSTRA, R. (1988). Order Restricted Statistical Inference. New York: Wiley.
TURNER, T. R. AND WOLLAN, P. C. (1997). Locating a maximum using isotonic regression. Computational Statistics and Data Analysis 25, 305–320.[CrossRef]
YEN, S. S. C. AND JAFFE, R. B. (1986). Reproductive Endocrinology: Physiology, Pathophysiology, and Clinical Management, 2nd edition. Philadelphia, PA: W. B. Saunders .
ZHANG, D., LIN, X., RAZ, J. AND SOWERS, M. (1998). Semiparametric stochastic mixed models for longitudinal data. Journal of the American Statistical Association 93, 710–719.[CrossRef][Web of Science]
ZHANG, D., LIN, X. AND SOWERS, M. (2000). Semiparametric regression for periodic longitudinal hormone data from multiple menstrual cycles. Biometrics 56, 31–39.[CrossRef][Web of Science][Medline]
Received February 11, 2004; revised January 26, 2005; accepted for publication February 2, 2005.
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||