PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Sankhya Ser B. Author manuscript; available in PMC 2010 October 16.
Published in final edited form as:
Sankhya Ser B. 2009; 71(1): 79–96.
PMCID: PMC2955899
NIHMSID: NIHMS202681

Statistical inference under order restrictions in analysis of covariance using a modified restricted maximum likelihood estimator

Abstract

In this article we introduce a new procedure for estimating population parameters under inequality constraints (known as order restrictions) when the unrestricted maximum liklelihood estimator (UMLE) is multivariate normally distributed with a known covariance matrix. Furthermore, a Dunnett-type test procedure along with the corresponding simultaneous confidence intervals are proposed for drawing inferences on elementary contrasts of population parameters under order restrictions. The proposed methodology is motivated by estimation and testing problems encountered in the analysis of covariance models. It is well-known that the restricted maximum likelihood estimator (RMLE) may perform poorly under certain conditions in terms of quadratic loss. For example, when the UMLE is distributed according to multivariate normal distribution with means satisfying simple tree order restriction and the dimension of the population mean vector is large. We investigate the performance of the proposed estimator analytically as well as using computer simulations and discover that the proposed method does not fail in the situations where RMLE fails. We illustrate the proposed methodology by re-analyzing a recently published rat uterotrophic bioassay data.

Keywords: Linear model, linked parameters, nodal parameter, order restrictions, restricted maximum likelihood estimators (RMLE), simple order, simple tree order, umbrella order

1 Introduction

In many applications researchers are interested in testing for inequality constraints among population means θi, i = 1, 2, …, p, after adjusting for covariates. For instance, toxicologists are often interested in studying the effect of a chemical on the mean weight of a specific organ of an animal after adjusting for its body weight (Kanno et al. 2002a, 2002b). Instead of the usual two-sided alternative Ha: θiθj, ij, researchers are often interested in testing against inequalities among the parameters (known as order restrictions). Some common order restrictions of interest (with at least one strict inequality) are; (a) Simple order θiθj, for all ij, (b) Simple tree order θ1θj, for all 1 ≤ j, (c) Umbrella order (with peak at i) θ1θ2 ≤ … ≤ θiθi+1 ≥ … ≥ θp. The null hypothesis being H0: θ1 = θ2 = … = θp.

Suppose [theta w/ hat]UMLE is the unrestricted maximum likelihood estimator (UMLE) of θ, with [theta w/ hat] ~ N(θ, σ2Σ), where Σ known. The problem of estimating θ under inequality constraints has received considerable attention in the literature. For an excellent review on this subject one may refer to the recent book by Silvapulle and Sen (2005). It is well known that the restricted maximum likelihood estimator (RMLE) of θ, under some inequality restrictions, may perform poorly in terms of mean squared error (MSE) and in terms of coverage probability for individual coordinates of θ (Lee, 1988 and Hwang & Peddada, 1994). Although the improved estimators derived in Hwang & Peddada (1994) perform well when Σ is diagonal, they may perform poorly for some non-diagonal patterns of Σ which do not satisfy the sufficient conditions described in their paper. The general iterative procedure developed in Peddada et al. (2005) can be equivalent or more efficient than the UMLE. The two procedures are equivalent when 1Σ−1 has some very small negative components, where 1 is a column vector of 1’s. For these reasons, in Section 2 we introduce a modification to RMLE which is demonstrated to perform better than the above estimators. It exploits the result that for p = 2 the RMLE universally dominates the UMLE (Rueda, et al., 1997). The proposed methodology is applicable for any partial order among the components of θ for p ≥ 2. Using this estimator, in Section 3 we introduce a Dunnett-type test procedure for testing the global null hypothesis of no difference among the components of θ against a partial order among θ. The proposed methodology is illustrated in Section 4 using the rat uterotrophic bioassay data discussed in Kanno et al. (2002b). Proofs are provided in the Appendix of the paper.

2 Estimation of parameters under order restrictions

2.1 Modified RMLE

Let Yij denote the response of the jth subject, j = 1, 2, …, ni, in the ith treatment group, i = 1, 2, …, p. Suppose xij denote a k × 1 vector of covariates associated with the jth subject in the ith treatment group. Then the standard analysis of covariance model is given by linear regression model

Yij=θi+xijβ+εij,i=1,2,,p,j=1,2,,ni,
(1)

where θi are the parameters of interest and β is the regression vector associated with the covariates xij and εij ~iid N(0, σ2). Rewriting the above equation in matrix notation, we have

YN×1=Diag(11,1p)θp×1+XN×kβk×1+εN×1,
(2)

where θ = (θ1, θ2, …, θp)′, 1i = (1, 1, …, 1)′ of order ni × 1, Diag is a block diagonal matrix with blocks made up of 1i, X is a N × k matrix with row vectors xij,N=i=1pni. Let [theta w/ hat]UMLE denote the UMLE of θ. Then, since εij ~iid N(0, σ2), we have that [theta w/ hat]UMLE ~ N(θ, σ2Σ), where Σ is a function of X and hence is known.

The basic idea of the proposed estimation procedure is the following: First, let us consider the special case when p = 2, and hence θ = (θ1, θ2)′. Assuming that [theta w/ hat]UMLE ~ N(θ, σ2Σ), the RMLE of θ is given by

θ^RMLE={θ^UMLEifθ^1UMLEθ^2UMLE(11θ^)(111)1(1,1)otherwise.
(3)

Recall that, this RMLE universally dominates the UMLE for this special case of p = 2 (see, Rudeda et. al. 1997). Our proposed method makes use of this result repeatedly. First, for each θi, reduce the p-dimensional order restricted problem to a collection of 2-dimensional problems in a way to described later. Then, for each of these 2-dimensional problems, we apply the foregoing special result for constructing an RMLE, and then average all such RMLEs of θi to obtain an ‘averaged’ estimator of θi (i = 1,, p). Since these averaged estimators may fail to satisfy the order restriction on θ, we use them as a set of initial estimates for constructing the final constrained estimator of θ to satisfy the order restrictions. The foregoing definitions make use of the graph theoretic ideas in Peddada et al (2005). Therefore, in what follows, we shall provide a brief description of the relevant graph theoretic ideas used in Peddada et al (2005) before describing the proposed method in detail.

A pair of parameters θi and θj is said to be linked if the inequality between them is known apriori. A parameter is said to be nodal if it is linked with every component of θ. There exists only one nodal parameter for simple tree and umbrella order restrictions, whereas every parameter is a nodal parameter in the case of simple order restriction.

For a parameter vector θ, the set of order restrictions is said to be a graph and is denoted by An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg. For a subvector of θ, the corresponding subset of order restrictions is said to be a subgraph. We shall use the phrases “(sub)graph”, “(sub)set” and “set of order restrictions” interchangeably. A subgraph x2133 is said to be a linked subgraph if all parameters in x2133 are linked, and it is said to be maximally linked subgraph if there does not exist a linked subgraph An external file that holds a picture, illustration, etc.
Object name is nihms202681ig2.jpg such that x2133 is a proper subset of An external file that holds a picture, illustration, etc.
Object name is nihms202681ig2.jpg. Thus in the case of simple order restriction θ1θ2 ≤ … ≤ θp, every subvector of θ corresponds to a linked subgraph and the maximally linked subgraph is An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg itself. Note that parameters belonging to a maximally linked subgraph satisfy a simple order restriction among the parameters belonging to the subgraph. If an order restriction satisfies inequalities of two maximally linked subgraphs x2133 and An external file that holds a picture, illustration, etc.
Object name is nihms202681ig2.jpg then we shall use the notation x2133ΛAn external file that holds a picture, illustration, etc.
Object name is nihms202681ig2.jpg to denote the order restriction. Suppose x21331, x21332, …, x2133m are the exhaustive collection of maximally linked subgraphs in An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg. Then C=Λi=1mMi. In the case of simple order restriction, we have m = 1 and An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg = x21331 = {θ [set membership] Rp: θ1θ2 ≤ … ≤ θp}, and for simple tree order we have m = p −1 with x2133i = {(θ1, θi) [set membership] R2: θ1θi}. In the case of umbrella order with peak at i, we have m = 2 with x21331 = {(θ1, θ2, …, θi) [set membership] Ri: θ1θ2 ≤ … ≤ θi−1θi} and x21332 = {(θi, θi+1, …, θp) [set membership] R(pi+1): θiθi+1 ≥ … ≥ θp}. We now define a lead parameter of a (sub)graph as the parameter that is linked to the largest number of parameters in the (sub)graph. If there is more than one such parameter then the lead parameter in the (sub)graph is the one whose UMLE has the smallest variance among the UMLEs of all such parameters. Thus in the case of simple order, the parameter whose UMLE has the smallest variance would be the lead parameter. In the case of umbrella order graph with peak at i, the lead parameter is θi and in the case of simple tree order graph it is θ1.

Let ω11, ω22, …, ωpp denote the diagonal elements of Σ, the covariance matrix of [theta w/ hat]UMLE. For a given component θi, let {θi1, θi2 …, θir} denote the set of all linked parameters (i = 1, 2, …, p). Using these UMLEs, for each linked pair (θi, θij)′, j = 1, 2, …, r, we obtain the RMLE using (3). Denote the RMLE by ( θ^i(i,j)RMLE,θ^ij(i,j)RMLE). Thus we obtain r different RMLE for θi, namely, θ^i(i,j)RMLE, j = 1, 2, …, r. Now, we average these RMLEs to obtain

θ^iave=j=1rαjθ^i(i,j)RMLEj=1rαj,

where αj is the inverse of ωijij, i.e. the inverse of the variance of θ^ijUMLE. Although each of the components used in the construction of θ^iave satisfy the desired order restriction, θ^iave, i = 1, 2, …, p may not satisfy the desired order restriction. To put them in the desired order, we begin with the lead parameter θk [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg and set θ^kθ^kave.

We now estimate all the parameters that are linked to θk. Identify all the maximally linked subgraphs which θk belongs to. Denote these by x2133k1, x2133k2, …, x2133kt. For each i = 1, 2, …, t, we estimate the parameters in x2133ki as follows. Suppose x2133ki is defined by θki1θki2 ≤ … θkil−1θkilθkil+1 ≤ … ≤ θkis, where θk [equivalent] θkil. Then the proposed order restricted estimators for θkij, j = 1, 2, …, s are:

θ^kij=min(θ^k(ij1)ave,θ^kij),j=l,l=1,,2andθ^kij=max(θ^kij,θ^kij+1ave),j=l,l+1,,s1.

We repeat the above process by identifying the next lead parameter in the complement of Λi=1tMki and repeat the above process until all the components of θ are estimated. Unlike the RMLE and the algorithm provided in Peddada et al. (2005), this procedure is not iterative.

Example 2.1

We now simplify the estimators for some common order restrictions to illustrate the above methodology.

  • Simple order (θ1θ2 ≤ … ≤ θp):
    Suppose the lead parameter is θk, where θ1θ2 ≤ … θk−1θkθk+1 … ≤ θp. Then θ^kθ^kave,θ^i=min(θ^iave,θ^i+1), i = k − 1, k − 2, …,1 and θ^i+1=max(θ^i,θ^i+1ave), i = k, k + 1, …, p − 1.
  • Simple tree order (θ1θi, i = 2, …, p):
    In this case the lead parameter is θ1. Therefore θ^1θ^1ave and for all i ≥ 2 θ^i=max(θ^1,θ^iRMLE).
  • Umbrella order (θ1θ2 ≤ … ≤ θk−1θkθk+1 ≥ … ≥ θp):
    In this case the lead parameter is θk. Therefore θ^kθ^kave,θ^i=min(θ^iave,θ^i+1) for i = k − 1, k − 2, …, 1, and θ^i+1=max(θ^i,θ^i+1ave) for i = k, k + 1, …, p − 1.

Theorem 2.1

Suppose [theta w/ hat]UMLE ~ N(θ, Σ) and suppose θ satisfies the order restrictions An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg. Then for all known positive constants αj, defined in θ^iave, and for every monotonic convex function ψ,

Eψ(θ^iaveθi)Eψ(θ^iUMLEθi),i=1,2,,p,
(4)

with a strict inequality holding for at least one i (i = 1, 2, …, p).

Relative to θ^iUMLE, if an estimator [theta w/ hat]i satisfies (4), then we shall say [theta w/ hat]i dominates θ^iUMLE in monotonic convex risk, or has smaller monotonic convex risk. Thus, by taking ψ(u) = u2 we deduce domination under the MSE criterion and by taking ψ (u) = |u| we deduce domination under the mean absolute error (MAE) criterion. Thus the domination result obtained in the above theorem is stronger than domination under the usual mean squared error and mean absolute error criteria.

As a consequence of the above theorem we have the following important corollary which justifies the proposed estimation procedure for any arbitrary order restriction An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg.

Corollary 2.2

Under the conditions of Theorem 2.1, suppose θk is the lead parameter of θ [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg then [theta w/ hat]k dominates θ^kUMLE in monotonic convex risk.

Remark 2.1

Under the conditions of Theorem 2.1, suppose θk is the only nodal parameter of θ [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms202681ig1.jpg then [theta w/ hat]k dominates θ^kUMLE in monotonic convex risk. Consequently we have the following important special cases:

  • For simple tree order θ1θi, i ≥ 2, [theta w/ hat]1 dominates θ^1UMLE in monotonic convex risk. Recall from Lee (1988) and Hwang and Peddada (1994) that the RMLE of θ1 fails to dominate the UMLE as p increases, whereas the proposed methodology yields an estimator which actually dominates the UMLE.
  • For umbrella order θ1θ2 ≤ … ≤ θk−1θkθk+1 ≥ … ≥ θp, [theta w/ hat]k dominates θ^kUMLE in monotonic convex risk.

2.2 Simulation studies

We evaluted the performance of the proposed modified RMLE (MRMLE) with the following estimators: (a) UMLE, (b) RMLE, (c) estimator of Hwang & Peddada (1994), and (d) the iterative algorithm of Peddada et al. (2005). In the simulation study reported in this paper, we assumed [theta w/ hat]~ N(θ, Σ), where the components of θ satisfy the simple tree order θ1θi, i = 2, 3, …, p, where θ1 is the mean of the control group and θi are the treatment means i ≥ 2. We considered the following patterns of Σ:

  • (P1) Independence: Σ = Diag[1: σ2: ···: σ2].
  • (P2) Intraclass among the p − 1 treatments: Σ = (σij, i, j = 1, 2, …, p), where σ11 = 1, σ1j = σ, j ≥ 2, σij = σ2, ij, i ≥ 2, j ≥ 2, σii = 1 + σ2, i ≥ 2.
  • (P3) Intraclass correlation: Σ = (1 − ρ)I + ρJ.

We considered two parameters, (1) the nodal parameter θ1, and (2) vector of elementary contrasts η = (η1, η2, …, ηp−1), where ηi = θi+1θ1. For any parameter η the MSE of an estimator [eta w/ hat] is defined as E([eta w/ hat]η)′([eta w/ hat]η). The coverage probability of a fixed width simultaneous confidence interval for η centered at [eta w/ hat] is defined as P(|[eta w/ hat]1η1| < cv1, |[eta w/ hat]2η2|< cv2, …, |[eta w/ hat]p−1ηp−1| < cvp−1), where vi is the standard error of the UMLE of ηi. The constant c is chosen so that P(η^1UMLEη1<cv1,η^2UMLEη2<cv2,,η^p1UMLEηp1<cvp1)=0.95, where η^iUMLE is the UMLE of ηi (i = 1, 2, …, p − 1). Thus for each estimator, the size of the region is the same, and hence, the estimator that has higher coverage probability is preferred as it has larger concentration of distribution around the true parameter. This is the peakedness criterion of Birnbaum (1948) for comparing estimators. We used these two criteria to evaluate the performance of an estimator. The univariate versions of these two criteria for the nodal parameter θ1 are defined accordingly.

All our simulation results are based on 10,000 simulation runs, with a variety of patterns for p, σ, ρ. Although a wide range of patterns were considered in our simulation study, in this paper we provide a representative sample of our results. Table 1 summarizes the results for the nodal parameter θ1 and Table 2 summarizes the results for the elementary contrast vector η. In all our simulations we considered the worst case scenario for θ, which is the boundary point (0, 0, …, 0).

Table 1
Mean Squared Error (MSE) and coverage rate comparisons of (a) UMLE, (b) RMLE, (c) estimator of Hwang & Peddada (1994), (d) the iterative algorithm of Peddada et al. (2005), and (e) the proposed modified RMLE (MRMLE). The parameter of interest ...
Table 2
Mean Squared Error (MSE) and coverage rate comparisons of (a) UMLE, (b) RMLE, (c) estimator of Hwang & Peddada (1994), (d) the iterative algorithm of Peddada et al. (2005), and (e) the proposed modified RMLE (MRMLE). The parameter of interest ...

As expected, in each case as the dimension increased, the RMLE performed very poorly in terms of MSE and the coverage probability. In several cases the MSE of RMLE was substantially larger than all its competitors. It also had poor coverage probability compared to the nominal level of 0.95. As expected, the estimator of Hwang & Peddada (1994) performed extremely well when Σ was diagonal (pattern P1) or when Σ had an intraclass correlation structure (pattern P3). However, its performance can be very poor for other covariance structures (pattern P2). The proposed estimator MRMLE consistently performed very well for all patterns considered in our simulation studies. Relative to RMLE, in some cases, its MSE was almost one-third or even one-fourth that of RMLE. Also, its coverage rate always exceeded the nominal level. Furthermore, as expected, the performance of MRMLE and the estimator proposed in Peddada et al. (2005) performed similarly for diagonal and intraclass correlation matrices. However, for other non-diagonal covariance matrices, the MRMLE performed better than the estimator in Peddada et al. (2005) in terms of MSE. In conclusion, MRMLE never failed the way RMLE and the estimator of Hwang & Peddada (1994) did, and was always better than the UMLE.

3 Hypothesis Testing

3.1 Dunnett-type procedure

The coverage probability results reported in Table 2 suggest that MRMLE may be used for constructing simultaneous confidence intervals for performing inferences on the elementary contrasts when the components of θ are subject to order restrictions. Analogous to Williams (1977) and Marcus and Talpaz (1992), we begin with the following testing problem along the lines of Dunnett (1955).

H0:θ1=θ2==θp,versusHa:θΘ,
(5)

where Θ denotes the set of order restrictions on the components of θ = (θ1, θ2, …, θp)′. Williams (1977) modified Dunnett (1955) when Θ is a simple order, whereas Marcus & Talpaz (1992) modified Dunnett’s procedure for simple tree order. In each case the basic idea is to replace the numerator of Dunnett’s test by the suitable RMLE. We adopt the same strategy, except that we shall use MRMLE in place of the RMLE as described below. This method generalizes the procedure provided in Peddada et al. (2006).

For an elementary contrast θiθj, the UMLE is θ^iUMLEθ^jUMLE. An estimator of the variance of this estimator is given by Var^(θ^iUMLEθ^jUMLE)=σ^2ai,jai,j, where ai,j is a vector of the form (0, 0, 1, …, −1, 0, …, 0)′ with +1 at the ith location and −1 at the jth location so that ai,jθ=θiθj. Let θ^i=fiΘ(θ^UMLE) be the MRMLE of θi. Then, analogous to Dunnett’s test, the proposed test statistic for testing (5) is given by

G=maxi,jθ^iθ^jσ^ai,jai,j,
(6)

where maxi,j is the maximum over all possible linked pairs i and j defined in Θ and [sigma with hat] is the usual mean sum of squares due to error.

Example 3.1

  1. Simple Order: If Θ is a simple order then (6) simplifies to
    G=max1ip1θ^i+1θ^iσ^ai+1,iai+1,i,
  2. Simple Tree Order: If Θ is the simple tree order restriction then (6) simplifies to
    G=max2iθ^iθ^1σ^ai,1ai,1.

Let Z ~ N(0, Σ) and let U ~ χ2 with df = Npk degrees of freedom. Under the order restriction specified in Ha, for each i, construct the MRMLE using Z. Denote it by Wi=fiΘ(Z). Let s2 denote the usual mean residual sum of squares in the linear model. Then we have the following theorem.

Theorem 3.1

Under the null hypothesis θ1 = θ2 = ··· = θp = θ and model (2)

maxi,jθ^iθ^jsai,jai,j=dmaxi,jdfWiWjUaijaij
(7)

In view of the above theorem, the null distribution of G can be approximated by simulating the following statistic:

G=maxi,jdfWiWjUai,jai,j.

Thus

P(GGαH0)=P(maxi,jθ^iθ^jσ^ai,jai,jGαH0)=α.

Hence we reject the global null H0 at a level of significance α if, for all linked parameters θi and θj,

θ^iθ^jσ^ai,jai,jGα.

Hence, for example if Ha is the simple tree order as defined above, then we reject H0: θ1 = θi, i ≥ 2, at a level of significance α if

maxi2θ^iθ^1σ^ai,1ai,1Gα.

As done in Dunnett (1955) and Marcus & Talpaz (1992) (see also Hsu, 1996), the above critical values can be used for constructing (1 − α) × 100% simultaneous confidence intervals for elementary contrasts θiθj between linked parameters θi and θj. Thus the simultaneous confidence intervals are given by

i,j{θiθjθ^iθ^j±Gασ^ai,jaij},θjθi.
(8)

The size of the above confidence hypercube is given by i,j{2×Gα×σ^ai,jai,j}.

3.2 Simulation Studies

We performed extensive simulation studies to compare the performance of the proposed procedure with the RMLE based procedure of Marcus & Talpaz (1992) when the parameters are subject to simple tree order restriction. We compared them in terms of (a) power and (b) size of the simultaneous confidence intervals for the p − 1 vector (θ2θ1, θ3θ1, …, θpθ1)′, where the size of the hypercube is defined earlier. For convenience, we report the natural logarithm of size. As in Section 2.2, we considered the patterns (P1) and (P3) for the covariance matrix of UMLE. The powers of the two procedures are summarized in Table 3 for the case (P1) since Marcus & Talpaz procedure was defined only for the diagonal covariance matrix. The size of simultaneous confidence intervals are provided in Table 4 for both covariance matrix patterns.

Table 3
Power comparison of Marcus & Talpaz test (MT) and the proposed test. Data are generated according Yij ~independent N (θi, 1), i = 1, 2, …, p, j = 1, 2, …, n, θ1θi. The powers are simulated using ...
Table 4
Size (log volumes) of simultaneous confidence intervals. The results are based on 10,000 simulation runs. The critical constants are chosen so that the coverage rate of the simultaneous confidence interval is 0.95.

From Table 3 it is apparent that the proposed procedure competes very well with MT procedure, for the extreme points such as θ = (0, 0, 0, …, 1)′, the MT procedure has higher power than the proposed test. However, for more central values such as θ = (0, 0, 1, 1, …, 1)′, the proposed method is more powerful than the RMLE based procedure of Marcus & Talpaz (1992).

Our simulation results reported in Table 4 suggest that the size of the simultaneous confidence intervals centered at MRMLE are smaller than those centered at RMLE. We observed this result for a variety of patterns. Hence the test based on the proposed simultaneous confidence intervals has larger power than the corresponding test based on RMLE.

4 Application to the study of estrogen-like compounds

Recently the Organization for Economic Cooperation and Development conducted a large uterotrophic bioassay to evaluate the effects of estrogen-like compounds on rat uterine weight (Kanno et al., 2002a, 2002b). One of the participating labs studied the effects of (1) Bisphenol A, (2) DBP, (3) DDT, (4) Genestein, (5) Ethinyl Estradiol (EE high dose), (6) Ethinyl Estradiol (EE low dose), (7) Methoxychlor.

Their study protocol consisted of 6 immature animals per test group, where the compound was administered using oral gavage for three consecutive days and the uterine weights were obtained on day 4 of the study. Immature animals were used so that the only source of estrogen was from the compounds used in the study.

The sample mean uterine weights and standard errors (in parentheses) along with mean body weight necropsy and the corresponding standard errors (in parentheses) are summarized in the following Table 5. Although the lab also studied the effects of Nonylphenol, we did not include this chemical in the present illustration since only 2 animals out of 6 survived in this group at the end of the study.

Table 5
Raw mean uterine weights and body weights at necropsy (day 4). Standard errors are in parentheses.

Since the uterine volume had a skewed distribution we log-transformed the data. Thus θi, i = 0, 1, …, 7, represents the mean uterine weight of log-transformed weights.

To illustrate the proposed methodology, we compute 95% simultaneous confidence intervals for θiθ0, i ≥ 1, under the hypothesized order restriction that θ0θi, i ≥ 1 with at least one strict inequality. If a particular confidence interval includes 0 then we conclude that the mean uterine weights of the particular treatment group is not significantly different from that of the control group.

Based on a two-way ANOVA, we found that the body weight was a significant covariate (P = 0.0238). Hence when comparing the above seven compounds with control group, we adjusted for the body weight of the animal. The unrestricted maximum likelihood estimates of mean uterine weights (in log scale) θi, i ≥ 0 (standard errors in parentheses) with body weight at necropsy as a covariate are as follows: 2.90 (0.12), 3.31 (0.12), 2.80 (0.10), 4.35 (0.12), 3.99 (0.11), 4.33 (0.10), 3.79 (0.10), 3.98 (0.10), respectively.

Under the simple tree order restriction, θ1θi, ∀i ≥ 1, the MRMLE are given by: [theta w/ hat]1 = 2.89, [theta w/ hat]2 = 3.31, [theta w/ hat]3 = 2.89, [theta w/ hat]4 = 4.35, [theta w/ hat]5 = 3.99, [theta w/ hat]6 = 4.33, [theta w/ hat]7 = 3.79, [theta w/ hat]8 = 3.98. Using the simulated distribution of G* we obtained a critical value of 2.42 for constructing the 95% simultaneous confidence intervals for θiθ1, i ≥ 2. Thus we have the following confidence intervals for θiθ1, i ≥ 2, 0.42 ± 2.42 × 0.17, 0 ± 2.42 × 0.17, 1.46±2.42×0.16, 1.1±2.42×0.16, 1.44±2.42×0.16, and 0.9±2.42×0.16, respectively.

Note that the 95% confidence interval for the mean change in uterine weights (relative to vehicle control) for the DBP group contains 0, suggesting no significant difference in the mean uterine weights between the DBP treated animals and the controls. None of the remaining confidence intervals contain 0. Thus the mean uterine volumes for all other test compounds differed from that of the control group. Except for Bisphenol A, our results regarding all other compounds agree with the results of Kanno et al. (2002b). They did not find Bisphenol A to be significantly different from the placebo group, but the proposed procedure did. Interestingly, our finding is in agreement with data obtained by some of the other labs that tested this compound (Kanno et al., 2002b).

Acknowledgments

The research of SDP was supported by the Intramural Research Program of the NIH, National Institute of Environmental Health Sciences[Z01 ES101744-04]. The authors thank Drs. Gregg Dinse and Anastasia Ivanova for carefully reading this manuscript and providing several useful comments that improved the presentation of the manuscript. We also thank the associate editor and the two referees for several insightful comments that led to a significant improvement in the presentation of the paper.

6 Appendix

Using the notations in (2) let

U=[Diag(11,1p):X]N×(p+k),

where the rank of U is p + k. Then the following lemma is useful for computing the covariance matrix of [theta w/ hat]UMLE.

Lemma 6.1

Suppose Ji is a ni × ni matrix of 1’s, M=Diag(J1n1,,Jpnp),L=Diag(1n1,,1np),1ni=(1/ni,1/ni,1/ni) is a ni × 1 vector, and E = (X′(IM)X)k×k, then

(UU)1=(Diag(1n1,,1np)+LE1XLLXE1E1XLE1)

Proof

Since Diag(n1, …, np) and XX are symmetric positive definite matrices and

UU=(Diag(n1,,np)Diag(11,,1p)XXDiag(11,1p)XX)

therefore from Rao (1973) we have

(UU)1=(Diag(1n1,,1np)+LE1XLLXE1E1XLE1)

Hence we have

Cov(θ^UMLE)=σ2=σ2(Diag(1n1,,1n1)+LXE1XL).

Accordingly, the variance of the UMLE of a contrast θiθj is given by

Var(θ^iUMLEθ^jUMLE)=σ2aijaij.
(9)

Proof of Theorem 2.1

From Rueda et al. (1997), we note that θ^i(i,j)RMLE universally dominates θ^iUMLE. Thus for all monotonic functions [var phi](·), E(φ(θ^i(i,j)RMLEθi))E(φ(θ^iRMLEθi)), with a strict inequality for at least one value of θi. Appealing to Dunbar et al. (2001) we therefore have E(ψ(θ^iaveθi))E(ψ(θ^iRMLEθi)), for all monotonic convex functions ψ(·). Hence the theorem.

Proof of Theorem 3.1

From (9) we have

maxi,jθ^iθ^jsVar(θ^iUMLEθ^jUMLE)=maxi,jσ(θ^iθ)(θ^jθ)saijaij.

Under the null hypothesis θ1 = θ2 = ··· = θp = θ we have θ^iθ^j=dσ(WiWj). Hence

maxi,jσ(θ^iθ)(θ^jθ)sσaijaij=dmaxi,jσWiWjsaijaij=dmaxi,jdfWiWjUaijaij

References

  • Birnbaum On random variables with comparable peakedness. Ann Math Statist. 1948;19:76–81.
  • Dunbar S, Conaway M, Peddada SD. On Improved Estimation of Parameters Subject to Order Restrictions. Statist and Appln. 2001;3:121–128.
  • Dunnett C. A multiple comparisons procedure for comparing several treatments with a control. J Amer Statist Assocn. 1955;50:1096–1121.
  • Fernandez M, Reuda C, Salvador B. The loss of efficiency estimating linear functions under order restrictions. Scand J Statist. 1999;26:579–592.
  • Fernandez M, Reuda C, Salvador B. Parameter estimation under orthant restriction. Canadian J Statist. 2000;28:171–181.
  • Hsu J. Multiple Comparisons: Theory and Methods. Chapman & Hall; New York, NY: 1996.
  • Hwang JTG, Peddada SD. Confidence Interval Estimation Subject to Order Restrictions. Ann Statist. 1994;22:67–93.
  • Kanno J, Onyon L, Peddada SD, Ashby J, Jacob E, Owens W. The OECD program to validate the uterotrophic bioassay: Phase Two - Dose Response Studies. Env Health Persp. 2002a;111:1530–1549. [PMC free article] [PubMed]
  • Kanno J, Onyon L, Peddada SD, Ashby J, Jacob E, Owens W. The OECD program to validate the rat uterotrophic bioassay: Phase Two - Coded Single Dose Studies. Env Health Persp. 2002b;111:1550–1558. [PMC free article] [PubMed]
  • Lee CIC. The Quadratic Loss of Order Restricted Estimators for Several Treatment Means and a Control Mean. Ann of Statist. 1988;16:751–758.
  • Marcus R, Talpaz H. Further results on testing homogeneity of normal means against simple tree alternatives. Comm Statist A. 1992;21:2135–2149.
  • Peddada SD, Dunson DB, Tan X. Estimation of order-restricted means from correlated data. Biometrika. 2005;92:703–715.
  • Peddada SD, Haseman J, Tan X, Travlos G. Tests for simple tree order restriction with application to dose-response studies. J Royal Statist Soc, Ser - C. 2006;55:493–506.
  • Rao CR. Linear Statistical Inference and Its Applications. Wiley; New York, NY: 1973.
  • Reuda C, Salvador B, Fernandez M. Simultaneous estimation in a restricted linear model. J Mult Anal. 1997;61:61–66.
  • Silvapulle M, Sen PK. Constrained Statistical Inference: inequality, order and shape restrictions. Wiley; New York, NY: 2005.
  • Williams D. Some inference procedures for monotonically ordered normal means. Biometrika. 1977;64:9–14.
  • Wollan P, Dykstra R. Algorithm AS 225: Minimizing Linear Inequality Constrained Mahalanobis Distances. Appl Stat. 1987;36:234–240.