Isotonic regression for metallic microstructure data: estimation and testing under order restrictions

ABSTRACT

Investigating the main determinants of the mechanical performance of metals is not a simple task. Already known physically inspired qualitative relations between 2D microstructure characteristics and 3D mechanical properties can act as the starting point of the investigation. Isotonic regression allows to take into account ordering relations and leads to more efficient and accurate results when the underlying assumptions actually hold. The main goal in this paper is to test order relations in a model inspired by a materials science application. The statistical estimation procedure is described considering three different scenarios according to the knowledge of the variances: known variance ratio, completely unknown variances, and variances under order restrictions. New likelihood ratio tests are developed in the last two cases. Both parametric and non-parametric bootstrap approaches are developed for finding the distribution of the test statistics under the null hypothesis. Finally an application on the relation between geometrically necessary dislocations and number of observed microstructure precipitations is shown.

Keywords:

• Isotonic regression
• alternating iterative method
• likelihood ratio test
• bootstrap
• order restrictions
• geometrically necessary dislocations

Abbreviations

AIM; alternating iterative method

EBSD; electron backscatter diffraction

GND; geometrically necessary dislocations

KAM; Kernel average misorientation

LRT; likelihood ratio test

MLE; maximum likelihood estimation

PAVA; pool-adjacentviolators algorithm

1. Introduction

Understanding the intrinsic nature of the mechanical properties of metals is usually not an easy task. In order to get insight into what gives desired mechanical performance to a metal, a deep and detailed analysis of the metal microstructure characteristics is needed. For instance, it is known in literature that dislocations, i.e. line defects in the crystalline arrangement of the atoms [17], play a fundamental role in the mechanical behavior of metal alloys. More specifically, the appearance of geometrically necessary dislocations¹ (GNDs) during plastic deformation of the material contributes to the hardening of the material. Detecting GNDs from 2D microstructure images is often challenging. One widely accepted way is to use the so-called Kernel average misorientation (KAM) [26]. The KAM, measured in electron backscatter diffraction (EBSD), quantifies the average misorientation around a measurement point with respect to a defined set of a nearest or nearest plus second-nearest neighbor points [8].

In [6,22,31], studies on the relation between GNDs and microstructure properties such as grain size and carbides size are presented. The relation between GNDs and grain size has both theoretical and experimental confirmation and it can be related to the well-known macroscopic physical Hall-Petch relation [15,29]. In fact, the Hall-Petch relation, in its original version, describes the negative dependence of yield stress (mechanical property) on grain size; loosely speaking the smaller grains are, the stronger the material is. More specifically in [18], the authors give as an explanation of the relation between GNDs and grain size that as the grain size decreases the grain boundary layer in which GNDs typically accumulate, occupies a greater volume fraction of the material, therefore it is reasonable to think that the smaller are the grains, the more GNDs will be observed.

Still unclear is, instead, the relation between carbides and GNDs. In fact, since the 1940s, several studies on how carbides affect the mechanical behavior of metals have been conducted. In [30], the authors state that the primary carbides and their distribution have a major influence on the wear resistance and the toughness of the material. However, carbides tend to precipitate along the grain boundaries, that as said before, are the locations in which GNDs typically accumulate. Until now, no direct physical relationship has been found between carbides and GNDs. Therefore, isolating carbides effect and assessing the conjecture on the positive relation between carbides and GNDs is a problem of interest.

In [16], a descriptive statistical analysis with response variable KAM, used as a proxy of GNDs and as explanatory variables the number of grains, the number of carbides and the position of carbides revealed an almost monotone trend of the response variable according to the increments of the explanatories.

Therefore, in order to take into account the already known direction of the physical relation, we want to propose an approach that incorporates this information and a procedure for testing the prementioned conjectures on a new dataset.

In this context, isotonic regression comes to aid. In fact, the idea at the basis of isotonic regression is taking order restrictions into account for improving the efficiency of the statistical analysis by reducing the error or the expected error of estimates and increasing the power of the testing procedures, provided that the hypothesized order restriction actually holds. The first papers about isotonic regression appeared in the 1950's [1,41] and books [2,33] are well-known references for frequentist statistical inference under order restrictions. Isotonic regression proves its power in different fields such as epidemiology in testing the effects of different treatments or in dose-finding [34,40] but also in genetics [23], business [19], biology [3]. There are not many examples of isotonic regression use in Materials Science. Throughout this paper, special attention is given to the peculiar data structure. Nowadays, developments toward multivariate isotonic regression, isotonic regression in inverse and censoring problems [13,14] are ongoing. Also in the Bayesian framework, the problem of estimating and testing under order constraints has been addressed and well developed [7,10,20,21]. But also in the most basic frequentist framework, there is still something missing.

In this paper, starting off with the most basic case, univariate isotonic regression of means under normality assumptions with known variances, we guide the reader into frequentist estimation and testing order restriction assumptions, considering different conditions on the variances.

Three different scenarios are considered. In all three cases, we focus on maximum likelihood as estimation procedure and likelihood ratio test (LRT) as test statistic for hypothesis testing.

The first case is the basic case in which ‘the variances’ are known or unknown but their ratio is known. This instance is considered extensively in [2,33] and results for estimation and testing order restrictions are already known.

The second scenario is from an applications point of view the most common scenario in which the variances are unknown. In [37], the authors derive a two steps estimating procedure for means and variances and interesting results on existence and uniqueness of the maximum likelihood estimates are derived under special conditions. Another iterative method, proposed in [38], is extended to the unknown variances case. The derivation of the test statistic and of its distribution in this scenario is not trivial. In fact, the estimate of the mean under the null hypothesis is also affected by the non-knowledge of the variances. We propose the LRT statistic and two different bootstrap approaches, one parametric and one non-parametric, for obtaining the test statistic distribution.

The last model considers not only the means under order restrictions but also the variances. This case has not often been faced probably because it is not common to have prior knowledge on the order of both means and variances. As in the unknown variances scenario, a two-step procedure for estimating means and variances is derived in [36] and similar results on existence and uniqueness under specific conditions on the empirical variances are given. In [38], an improved algorithm called alternating iterative method (AIM) and more general results about convergence are derived. For testing in this case we derive the LRT taking into account the order of variances also under the null hypothesis and apply a parametric and non-parametric bootstrap approach in line with the one derived in the unknown variance case to obtain approximate p-values.

The paper structure is as follows. In Section 2, we explain the estimation procedure of the isotonic means in the three different cases (Sections 2.1, 2.2, 2.3). In Section 3, the focus is on the LRT. We present it in the three different cases (Sections 3.1, 3.2, 3.3) and in Section 4, we propose both a parametric and non-parametric bootstrap approach for approximating the distribution of the test statistics under the null hypothesis. Finally, in Section 5, we come back to the application and we illustrate step-by-step how to deal with a real problem and more precisely how to perform isotonic regression and test for monotonicity of KAM with respect to the number of carbides. The paper ends with conclusions in Section 6.

2. Estimating restricted means in the normal case

We first introduce isotonic regression and the notation used in the rest of the paper in a more general context. Normality is assumed throughout this section.

Let $y_{ij}$, $j=1,\dots ,n_i$, $i=1,\dots ,k$ be the jth observation of the response variable $\mathrm{Y}$ corresponding to the ith level of the explanatory variable $\mathrm{X}$.

We assume ${\mathrm{Y}}_{ij}$ to be independent random variables, normally distributed with means ${\mu }_i$ and variances ${\sigma }_i^2$, $i=1,\dots ,k$, $j=1,2,\dots ,n_i$.

The log-likelihood is then given by (1) $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2\right\}+\mathrm{c}$(1) where $\mathrm{c}$ is a constant which does not depend on the parameters $\mathit{\mu }=({\mu }_1,\dots {\mu }_k{)}^{\prime }$ and ${\mathit{\sigma }}^{\mathbf{2}}=({\sigma }_1^2,\dots ,{\sigma }_k^2{)}^{\prime }$.

Furthermore, we assume that $\mathit{\mu }$ satisfies (2) ${\mu }_1\le {\mu }_2\le \cdots \le {\mu }_k.$(2) A k-dimensional vector $\mathit{\mu }$ is said to be isotonic if $t\le s$ implies ${\mu }_t\le {\mu }_s$.

Let D be the set of all the isotonic vectors in ${\mathbb{R}}^k$, (3) $D=\{\mathit{\mu }\in {\mathbb{R}}^k;{\mu }_1\le {\mu }_2\le \cdots \le {\mu }_k\}$(3) In this section, we are interested in the maximum likelihood estimator of $(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$, where $\mathit{\mu }$ is isotonic and ${\sigma }_i^2>0$. Depending on the information on ${\mathit{\sigma }}^{\mathbf{2}}$, different MLEs have been derived.

In the following three sections, the three different cases are considered.

2.1. Isotonic regression of means with known variance ratio

This first case constitutes the most basic case in which all variances are either known or unknown but they differ according to some known multiplicative constants $c_i$. This means that the variance ${\sigma }_i^2$ of the response variable ${\mathrm{Y}}_i$ is given by: ${\sigma }_i^2=c_i{\sigma }^2,1\le i\le k.$ This specific case is already covered in [2,33], but we hereafter report the main results. The problem of maximizing log-likelihood (1(1) $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2\right\}+\mathrm{c}$(1) ) in μ can be rewritten equivalently as solving: (4) $\underset{\mu \in D}{min}\sum _{i=1}^k({\overline{y}}_i-{\mu }_i{)}^2{w}_i$(4) where ${\overline{y}}_i=\frac{\sum _j{y}_{ij}}{n_i}$ and $w_i=\frac{n_i}{c_i}$. Note that this objective function does not depend on ${\sigma }^2$. The solution, ${\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}^I$, is called the isotonic regression of $\overline{\mathit{y}}=({\overline{y}}_1,\dots ,{\overline{y}}_k)$ with weights $\mathit{w}=(w_1,\dots ,w_k)$ [37]. For obtaining the solution to (4(4) $\underset{\mu \in D}{min}\sum _{i=1}^k({\overline{y}}_i-{\mu }_i{)}^2{w}_i$(4) ), different algorithms have been proposed in the literature ([2,33]). In this paper, the the intuitive and skilfully implemented ‘Pool-Adjacent Violators Algorithm’ (PAVA) is used ([1,5]).

More details about the algorithm are provided in Appendix 1.

2.2. Isotonic regression of means with unknown variances

In this second case, no assumptions on the variances are made. They are unknown and for obtaining the maximum likelihood estimate of $\mathit{\mu }$, they need to be estimated as well. In [37], the authors consider this case and interesting results on existence and uniqueness of the MLE are achieved. We hereby recall the main results. The approach is to maximize the log-likelihood (1(1) $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2\right\}+\mathrm{c}$(1) ), with $\mathit{\mu }\in D$ and ${\mathit{\sigma }}^{\mathbf{2}}\in {\mathbb{R}}_{+}^k$.

For any fixed ${\mathit{\sigma }}^{\mathbf{2}}\in {\mathbb{R}}_{+}^k$, the maximizer ${\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}^{\mathit{I}}$ of $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$ over $\mathit{\mu }\in D$ is the isotonic regression of $\overline{\mathit{y}}$ with weights $\mathit{w}=(w_1,\dots w_k{)}^{\prime }$ and $w_i=\frac{n_i}{{\sigma }_i^2}$.

On the other hand, for any fixed $\mathit{\mu }\in D$, the maximizer ${\mathit{\sigma }}^{\mathbf{2}}$ of $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$ over $\mathit{\sigma }\in {\mathbb{R}}_{+}^k$ is ${\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}^{\mathbf{2}}\left(\mathit{\mu }\right)=\left({\hat{\sigma }}_1^2\right({\mu }_1),\dots {\hat{\sigma }}_k^2({\mu }_k){)}^{\prime }$, where ${\hat{\sigma }}_i^2\left({\mu }_i\right)=\frac{\sum _j^{n_i}(y_{ij}-{\mu }_i{)}^2}{n_i}$.

Substituting ${\hat{\sigma }}^2\left(\mu \right)$ into (1(1) $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2\right\}+\mathrm{c}$(1) ), we can express the profile log-likelihood of $\mathit{\mu }$ as (5) $l\left(\mathit{\mu }\right)=\sum _{i=1}^k-n_i\mathrm{l}\mathrm{n}[{\overline{\sigma }}_i^2+({\overline{y}}_i-{\mu }_i{)}^2]+c$(5) where ${\overline{\sigma }}_i^2=\frac{\sum _j^{n_i}(y_{ij}-{\overline{y}}_i{)}^2}{n_i}$ is the sample variance of the ith normal population and c a constant that does not depend on $\mathit{\mu }$. Note that $l\left(\mathit{\mu }\right)\to -\mathrm{\infty }$ if ${\mu }_k\to \mathrm{\infty }$ or ${\mu }_1\to -\mathrm{\infty }$. Hence, maximizing l over D is equivalent to maximizing l over a compact subset of D of type $D_a=\{\mu \in D:{\mu }_1\ge -a,{\mu }_k\le a\}$. As l is continuous on $D_a$, a maximizer over D exists.

As previously said, the authors in [37] discuss also uniqueness of the MLE of $(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$. They state that l is not a concave function in general and that for guaranteeing uniqueness the following condition suffices (see Theorem 2.3 [37]):

Condition 2.1

For $i=1,\dots ,k$, ${\overline{\sigma }}_i^2>max\left\{\right({\overline{y}}_i-min\left(\overline{\mathit{y}}\right){)}^2,({\overline{y}}_i-max(\overline{\mathit{y}}\left){)}^2\right\}$.

For finding a maximizer of (5(5) $l\left(\mathit{\mu }\right)=\sum _{i=1}^k-n_i\mathrm{l}\mathrm{n}[{\overline{\sigma }}_i^2+({\overline{y}}_i-{\mu }_i{)}^2]+c$(5) ), a two-step iterative algorithm based on PAVA has been proposed in [37]. From an initial guess for $\mathit{\mu }$, the associated maximizer in ${\mathit{\sigma }}^{\mathbf{2}}$ is computed and after that the maximizer in $\mathit{\mu }$ based on this ${\mathit{\sigma }}^{\mathbf{2}}$ and so on. This iterative procedure stops when the maximum difference between the estimated means at step l−1 and at step l is less than an arbitrary small threshold value, e.g. $\underset{1\le i\le k}{max}|{\mu }_i^{I(l-1)}-{\mu }_i^{I\left(l\right)}|\le {10}^{-m},$ where m is taken to be equal to 3 in our case. In [38], the authors propose a new algorithm called AIM. The procedure is based on the minimization of a semi-convex function. In particular, restating the problem in terms of $(\mathit{\mu },\mathit{\nu })$, where $\mathit{\nu }=(1/{\sigma }_1^2,\dots ,1/{\sigma }_k^2{)}^{\prime }$ and given $D_a$ is a convex subset of ${\mathbb{R}}^k$ and V a convex subset of ${\mathbb{R}}_{+}^k$, $V=\{\nu \in {\mathbb{R}}_{+}^k:0\le 1/\underset{i}{max}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right))\le {\nu }_i\le 1/\underset{i}{min}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)\left)\right\}$, $L(\mathit{\mu },\mathit{\nu })$ is a semi-convex function because: i) $L(\mathit{\mu },\mathit{\nu })$ is defined on $D_a\times V$; ii) for any given $\mathit{\mu }\in D_a$, $L(\mathit{\mu },\cdot )$ is strictly convex on V and, for any given $\mathit{\nu }\in V$, $L(\cdot ,\mathit{\nu })$ is strictly convex on $D_a$. The algorithm originally proposed for the simultaneous order restrictions of means and variances can be easily extended to the unknown variance case. The iteration method works in alternating the search of the minimum point, ${\mathit{\mu }}^{\left(l\right)}$, of $L(\mathit{\mu },\mathit{\nu }({\mathit{\mu }}^{(l-1)}\left)\right)$ on a compact subset $D_a$ and the search of the minimum point, ${\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$, of $L\left(\mathit{\mu }\right({\mathit{\nu }}^{(l-1)}),\mathit{\nu })$ on V. Proof of the convergence of the algorithm does not require additional conditions [38]. The iterative procedure stops when the difference between the likelihoods at step l−1 and at step l is less than an arbitrary small threshold value: (6) $|L({\mathit{\mu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}},{\mathit{\nu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}})-L(;{\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}},{\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}})|\le {10}^{-m}$(6) A more detailed version of both algorithms is reported in Appendix A.2.

2.3. Isotonic regression of means and variances simultaneously

We now assume that both mean and variances are restricted by simple orderings. Therefore, in addition to assumption (2(2) ${\mu }_1\le {\mu }_2\le \cdots \le {\mu }_k.$(2) ), we assume also: (7) ${\sigma }_1^2\ge {\sigma }_2^2\ge \cdots \ge {\sigma }_k^2>0$(7) The reason for taking decreasing order relates to our application considered in Section 5; increasing variances can be dealt with analogously. In [36], maximum likelihood estimation under simultaneous order restrictions on mean and variances from a normal population is studied. Some of the most important results are hereby recalled. The approach is to maximize the log-likelihood (1(1) $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2\right\}+\mathrm{c}$(1) ) with $\mathit{\mu }\in D$ and ${\mathit{\sigma }}^{\mathbf{2}}\in \overline{G}$, where $\overline{G}$ is the closure of (8) $G=\{{\mathit{\sigma }}^{\mathbf{2}}\in {\mathbb{R}}_{+}^k:{\sigma }_1^2\ge {\sigma }_2^2\ge \cdots \ge {\sigma }_k^2>0\}.$(8) This means that the maximizer will have positive ${\sigma }^2$-values if there is variation within the groups. Then, for any fixed ${\mathit{\sigma }}^{\mathbf{2}}\in G$, the maximizer ${\mathit{\mu }}^{\mathit{I}}$ of $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$ over $\mathit{\mu }\in D$ is the isotonic regression of $\overline{\mathit{y}}$ with weights $\mathit{w}=(w_1,\dots w_k{)}^{\prime }$ and $w_i=\frac{n_i}{{\sigma }_i^2}$.

Furthermore, for any $\mathit{\mu }\in D$, the maximizer ${\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}^{\mathbf{2}\mathit{I}}\left(\mathit{\mu }\right)$ of $l(\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})$ is the so-called antitonic regression (isotonic regression with reversed order [13]) of ${\mathit{s}}^{\mathbf{2}}=(s_1^2,\dots s_k^2{)}^{\prime }$, $s_i^2=\frac{\sum _{j=1}^{n_i}(y_{ij}-{\mu }_i{)}^2}{n_i}$, with weights $\mathit{N}=(n_1,\dots ,n_k{)}^{\prime }$. Existence is guaranteed noticing that ${\mathit{\sigma }}^{\mathbf{2}}\in \left[\underset{i}{min}\right(\underset{\underset{i}{min}\left(\overline{\mathit{y}}\right)\le \theta \le \underset{i}{max}\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)),\underset{i}{max}(\underset{\underset{i}{min}\left(\overline{\mathit{y}}\right)\le \theta \le \underset{i}{max}\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)\left)\right]$, $s_i^2\left(\theta \right)=\sum _{j=1}^{n_i}(y_{ij}-\theta {)}^2/n_i$ (see [36, Theorem 2.1]).

Uniqueness is proved under the following condition (see Theorem 2.2 [36])

Condition 2.2

For $i=1,\dots k$ the sample variance ${\overline{\sigma }}_i^2$ satisfies ${\overline{\sigma }}_i^2>2(b-a)$, where b and a are the maximal and the minimal means respectively.

As in the unknown variances case, a two-step iterative algorithm is proposed for finding the solution for both means and variances under order restrictions. The proof of the convergence of the algorithm is given under Condition 2.2.

Later, in [38], as mentioned in the previous section, the authors show that restating the problem in terms of $(\mathit{\mu },\mathit{\nu })$, where $\mathit{\nu }=(1/{\sigma }_1^2,\dots ,1/{\sigma }_k^2{)}^{\prime }$ Condition 2.2 is not needed for proving that the algorithm converges. In fact, also in this case the proposed AIM algorithm can be employed. Since $L(\mathit{\mu },\mathit{\nu })$ has continuous second-order partial derivatives and the Hessian matrix with respect to μ $H(\mathit{\mu },\mathit{\nu })=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(n_1{\nu }_1,\dots ,n_k{\nu }_k)$ is a positive definite diagonal matrix for any fixed $\mathit{\nu }=({\nu }_1,\dots ,{\nu }_k{)}^{\prime }\in V_0$, $V_0=\{\nu \in {\mathbb{R}}^k:0\le 1/\underset{i}{max}(\underset{\underset{i}{min}\left(\overline{\mathit{y}}\right)\le \theta \le \underset{i}{max}\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right))\le {\nu }_1\le \cdots \le {\nu }_k\le 1/\underset{i}{min}(\underset{\underset{i}{min}\left(\overline{\mathit{y}}\right)\le \theta \le \underset{i}{max}\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)\left)\right\}$ then by Theorem 4 in [38] the iterative sequence of solutions to $L(\mathit{\mu },\mathit{\nu })$, $\left\{\right({\mathit{\mu }}^{\mathbf{\left(}\mathit{n}\mathbf{\right)}},{\mathit{\nu }}^{\mathbf{\left(}\mathit{n}\mathbf{\right)}}\left)\right\}$ converges to the MLE solution and consequently the sequence $\left\{\right({\mathit{\mu }}^{\mathbf{\left(}\mathit{n}\mathbf{\right)}},{\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathit{n}\mathbf{\right)}}\left)\right\}$ as well.

As in the previous case, the alternating iterative procedure is stopped when the maximum difference between the likelihoods at step l−1 and at step l is less than an arbitrary small threshold value (see (6(6) $|L({\mathit{\mu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}},{\mathit{\nu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}})-L(;{\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}},{\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}})|\le {10}^{-m}$(6) )).

A pseudo code of the algorithms can be found in Appendix A.3.

3. LRT: constant μ against monotonicity

We are interested in testing hypotheses of monotonicity in μ under the various assumptions on the variances discussed in Section 2. There exists extensive literature on testing hypotheses on means. In most cases, a standard testing procedure entails testing the hypothesis of equality of means against the hypothesis that they are different. In this paper, we consider the same null hypothesis but the alternative is different: monotonicity of the means. As in the previous section, we consider three different testing frameworks according to the different assumptions on the variances. In all three different scenario, the test statistic of interest is the LRT, an intuitive and powerful tool in hypothesis testing. In both [2] and [33], an entire chapter is dedicated to LRT developments and its use for testing order restrictions hypothesis under the normality assumption and known variance ratio. Using the same notation used in Section 2, we wish to test $H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_k$ against monotonicity of means (9) $H_1:{\mu }_1\le {\mu }_2\le \cdots \le {\mu }_k.$(9) The LRT for $H_0$ against $H_1$ can be defined as: (10) $\mathrm{\Lambda }=\frac{\underset{(\mathit{\mu }\in H_0;{\mathit{\sigma }}^{\mathbf{2}})}{max}L({\mathit{y}}_{\mathbf{1}},{\mathit{y}}_{\mathbf{2}},\dots ,{\mathit{y}}_{\mathit{k}};\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})}{\underset{(\mathit{\mu }\in H_1;{\mathit{\sigma }}^{\mathbf{2}})}{max}L({\mathit{y}}_{\mathbf{1}},{\mathit{y}}_{\mathbf{2}},\dots ,{\mathit{y}}_{\mathit{k}};\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})}$(10) where ${\mathit{y}}_{\mathit{i}}=(y_{i1},\dots ,y_{i{n}_i}{)}^{\prime }$, $\mathit{\mu }=({\mu }_1,\dots {\mu }_k{)}^{\prime }$ and ${\mathit{\sigma }}^{\mathbf{2}}=({\sigma }_1^2,\dots {\sigma }_k^2{)}^{\prime }$. It rejects the null hypothesis for small values of Λ or alternatively for large values of $-2\log \mathrm{\Lambda }$. The convenience in using this other form lies on the analogy with the ${\chi }^2$ statistic used to test against the alternative hypothesis ${\overline{H}}_0$, that not all ${\mu }_i$'s, $i=1,\dots ,k$, are the same.

In the following sections, more explicit expressions for Λ are given depending on the specific assumptions on means and variances.

3.1. LRT with known variance ratio

As in Section 2.1, let $y_{ij}$ $j=1,2,\dots n_i$, $i=1,2,\dots k$ be independent observations, normally distributed with unknown mean ${\mu }_i$ and variances ${\sigma }_i^2=c_i{\sigma }^2$ with $c_i$ known and ${\sigma }^2$ unknown. Under $H_0$, the maximum likelihood estimate of ${\mu }_1={\mu }_2=\cdots ={\mu }_k$ is given by: (11) ${\hat{\mu }}_{H_0}=\frac{\sum _{i=1}^k{w}_i{\overline{y}}_i}{\sum _{i=1}^k{w}_i}$(11) with $w_i=\frac{n_i}{c_i}$. Under $H_1$ the MLE of $\mathit{\mu }$ is ${\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}_{{\mathit{H}}_{\mathbf{1}}}^{\mathit{I}}$, the isotonic regression of $\overline{\mathit{y}}$, with weights $\mathit{w}=(w_1,\dots ,w_k{)}^{\prime }$, with respect to the simple order defined in (9(9) $H_1:{\mu }_1\le {\mu }_2\le \cdots \le {\mu }_k.$(9) ).

The LRT for $H_0$ against $H_1$, if the variances are known and $c_i=1$ boils down to rejecting $H_0$ for large values of (12) $-2\log \mathrm{\Lambda }=\frac{1}{{\sigma }^2}\left[\sum _{i=1}^k\sum _{j=1}^{n_i}(y_{ij}-{\hat{\mu }}_{H_0}{)}^2-\sum _{i=1}^k\sum _{j=1}^{n_i}(y_{ij}-{\hat{\mu }}_{i{H}_1}^I{)}^2\right]$(12) It is easy to check that the test is equivalent to rejecting $H_0$ for large values of: (13) ${\overline{\chi }}^2=\frac{\sum _{i=1}^k{\overline{\chi }}_i^2}{{\sigma }^2}$(13) where ${\overline{\chi }}_i^2=n_i({\hat{\mu }}_{i{H}_1}^I-{\hat{\mu }}_{H_0}{)}^2$ and ${\sigma }^2$ is the (known) common value of the variance.

Now, let us consider the more general case, ${\sigma }_i^2=c_i{\sigma }^2$ with $c_1,c_2,\dots c_k$ known and ${\sigma }^2$ unknown. The estimator of ${\sigma }^2$ under the null hypothesis is (14) ${\hat{\sigma }}_{H_0}^2=\frac{\sum _{i=1}^k{c}_i^{-1}\sum _{j=1}^{n_i}(y_{ij}-{\hat{\mu }}_{H_0}{)}^2}{N}$(14) and under $H_1$ (15) ${\hat{\sigma }}_{H_1}^2=\frac{\sum _{i=1}^k{c}_i^{-1}\sum _{j=1}^{n_i}(y_{ij}-{\hat{\mu }}_{i{H}_1}^I{)}^2}{N}$(15) The LRT rejects $H_0$ for small values of $\mathrm{\Lambda }=(\frac{{\hat{\sigma }}_{H_1}^2}{{\hat{\sigma }}_{H_0}^2}{)}^{N/2}$ or equivalently, taking ${\overline{E}}^2=1-{\mathrm{\Lambda }}^{2/N}$, for large values of (16) ${\overline{E}}^2=\frac{\sum _{i=1}^k{c}_i^{-1}{\overline{\chi }}_i^2}{\sum _{i=1}^k{c}_i^{-1}\sum _{j=1}^{n_i}(y_{ij}-{\hat{\mu }}_{H_0}{)}^2}$(16) An extension to the multivariate case with covariance matrix Σ unknown but common can be found in [28,35].

3.2. LRT with unknown variances

In this second case, no assumptions on the variances are made. They are unknown and possibly unequal. Using the notation of Section 2.2, let $y_{ij}$, $j=1,2,\dots ,n_i$, $i=1,2,\dots ,k$ be independent observations from a univariate normal distribution with unknown mean vector ${\mu }_i$ and completely unknown variances ${\sigma }_i^2>0$. Let ${\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}^{\mathit{I}}$ be the solution of the isotonic regression of $\overline{\mathit{y}}$ with weights $\mathit{w}=(w_1,\dots ,w_k{)}^{\prime }$, $w_i=\frac{n_i}{{\sigma }_i^2}$ found used Algorithm (2.2) in Appendix 1.

The first example of testing when all the variances are unknown can be found in [4] and the univariate version of the test proposed by the author is: (17) $\sum _{i=1}^k\frac{({\hat{\mu }}_i^I-\overline{y}{)}^2{n}_i}{s_i^2}$(17) where $\overline{y}=\frac{\sum _{i=1}^k{n}_i{\overline{y}}_i}{\sum _{i=1}^k{n}_i}$ and $s_i^2=\frac{\sum _{i=1}^k\sum _{j=1}^{n_i}(y_{ij}-{\overline{y}}_i{)}^2}{n_i-1}$. This test is clearly inspired by the LRT but it is not.

Let us consider first the maximum likelihood solution $({\hat{\mu }}_{H_0},{\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}_{H_0}^{\mathbf{2}})$, ${\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}_{{\mathit{H}}_{\mathbf{0}}}^{\mathbf{2}}=({\hat{\sigma }}_{1H_0}^2,\dots ,{\hat{\sigma }}_{k{H}_0}^2{)}^{\prime }$ under the null hypothesis. The log-likelihood under the null hypothesis is (18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) Differentiating this log-likelihood with respect to μ and ${\sigma }_i^2$, the following k + 1 score equations in k + 1 unknowns emerge: (19) $\left\{\begin{array}{ll}{\mu }_{H_0}=\frac{\sum _{i=1}^k{n}_i{\sigma }_{i{H}_0}^{-2}{\overline{y}}_i}{\sum _{i=1}^k{n}_i{\sigma }_{i{H}_0}^{-2}}& \\ {\sigma }_{i{H}_0}^2=\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-{\mu }_{H_0}{)}^2& 1\le i\le k\end{array}\right.$(19) Substituting ${\mathit{\sigma }}_{{\mathit{H}}_{\mathbf{0}}}^{\mathbf{2}}\mathbf{\left(}\mathit{\mu }\mathbf{\right)}$ in (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ), the profile likelihood of μ is: (20) $l\left(\mu \right)=-\sum _{i=1}^k\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-\mu {)}^2\right)+c.$(20)

Theorem 1

A maximizer of (20(20) $l\left(\mu \right)=-\sum _{i=1}^k\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-\mu {)}^2\right)+c.$(20) ) over ${\mathbb{R}}^d$ exists and it is contained in $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$. Moreover, if $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]\in \left[\underset{1\le i\le k}{max}\right({\overline{y}}_i-{\overline{\sigma }}_i),\underset{1\le i\le k}{min}({\overline{y}}_i+{\overline{\sigma }}_i\left)\right]$ then the maximizer is unique.

Proof.

Maximizing profile likelihood of μ (20(20) $l\left(\mu \right)=-\sum _{i=1}^k\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-\mu {)}^2\right)+c.$(20) ) boils down to maximize the sum of functions (21) $-\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(n_i\right({\overline{\sigma }}_i^2+({\overline{y}}_i-\mu {)}^2)),i=1,\dots ,k.$(21) Functions of type (21(21) $-\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(n_i\right({\overline{\sigma }}_i^2+({\overline{y}}_i-\mu {)}^2)),i=1,\dots ,k.$(21) ) are unimodal with mode at ${\overline{y}}_i$ and strictly concave on $[{\overline{y}}_i-{\overline{\sigma }}_i;{\overline{y}}_i+{\overline{\sigma }}_i]$. As the sum of unimodal functions is decreasing to the right of the rightmost mode (since all terms are decreasing) and from $-\mathrm{\infty }$ to the leftmost mode, the sum is increasing (as all of the functions are increasing on that set). Therefore, any maximizer of l, if it exists, belongs to the interval $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$. As l is continuous on $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$, existence of a maximizer is guaranteed.

Then if we consider the (possibly empty) interval where all the functions in (21(21) $-\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(n_i\right({\overline{\sigma }}_i^2+({\overline{y}}_i-\mu {)}^2)),i=1,\dots ,k.$(21) ) are strictly concave, on that interval the sum is also strictly concave. As for each i the function (21(21) $-\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(n_i\right({\overline{\sigma }}_i^2+({\overline{y}}_i-\mu {)}^2)),i=1,\dots ,k.$(21) ) is strictly concave on $I_i=[{\overline{y}}_i-{\overline{\sigma }}_i;{\overline{y}}_i-{\overline{\sigma }}_i]$, (20(20) $l\left(\mu \right)=-\sum _{i=1}^k\frac{n_i}{2}\mathrm{l}\mathrm{n}\left(\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-\mu {)}^2\right)+c.$(20) ) is strictly concave on $\bigcap _{i=1}^k{I}_i$. If $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$ is contained in this intersection, l is strictly concave on $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$. Hence l has a unique maximizer on ${\mathbb{R}}^d$.

Remark 3.1

Remark in a setting with real data, it is easy to check whether $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]\in \left[\underset{1\le i\le k}{max}\right({\overline{y}}_i-{\overline{\sigma }}_i),\underset{1\le i\le k}{min}({\overline{y}}_i+{\overline{\sigma }}_i\left)\right]$ and hence to determine whe-ther the maximum is unique.

However, as seen from (19(19) $\left\{\begin{array}{ll}{\mu }_{H_0}=\frac{\sum _{i=1}^k{n}_i{\sigma }_{i{H}_0}^{-2}{\overline{y}}_i}{\sum _{i=1}^k{n}_i{\sigma }_{i{H}_0}^{-2}}& \\ {\sigma }_{i{H}_0}^2=\sum _{j=1}^{n_i}{n}_i^{-1}(y_{ij}-{\mu }_{H_0}{)}^2& 1\le i\le k\end{array}\right.$(19) ), the MLE estimate $(\mu ,{\mathit{\sigma }}^{\mathbf{2}})$ has no closed form expression. Therefore, in [11] and [27], two different methods for finding the optimal solution are proposed. The first is an iterative procedure based on the Newton-Raphson method. A reasonable initial value for ${\hat{\mu }}_{H_0}^{\left(0\right)}$ is the so-called Graybill-Deal estimator [12] ${\hat{\mu }}_{\left(GD\right)}=\frac{\sum _{i=1}^k\left(n_i{\overline{y}}_i\right)/{\overline{s}}_i^2}{\sum _{i=1}^k{n}_i/{\overline{s}}_i^2}$ with ${\overline{s}}_i^2=\frac{\sum _{j=1}^{n_i}(y_{ij}-{\overline{y}}_i{)}^2}{n_i-1}$. The convergence speed of the algorithm strongly depends on the initial values. The second method is based on the profile likelihood approach. The authors in [27] propose the bisection method for finding the zero of the profile likelihood with respect to ${\mu }_{H_0}$. Under $H_1$ we use as estimates of $({\mathit{\mu }}_{{\mathit{H}}_{\mathbf{1}}},{\mathit{\sigma }}_{\mathit{i}\mathit{H}\mathbf{1}}^{\mathbf{2}})$, $({\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}^{\mathit{I}},{\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}^{\mathbf{2}})$ found using the iterative procedure described in Section 2.2.

The LRT when the variances are completely unknown can be expressed as: $\tilde{\mathrm{\Lambda }}=\prod _{i=1}^k{\left(\frac{{\hat{\sigma }}_{i{H}_0}^2}{{\hat{\sigma }}_{i{H}_1}^2}\right)}^{-\frac{n_i}{2}}$ Therefore, as in the previous case, the test rejects for small values of $\tilde{\mathrm{\Lambda }}$ or equivalently for large values of $-2\log \tilde{\mathrm{\Lambda }}$.

3.3. LRT with ordered variances

Using the notation of Section 2.3, let $y_{ij}$, $j=1,2,\dots ,n_i$, $i=1,2,\dots ,k$ be independent observations from Normal distributions with mean vector ${\mu }_i$ and variances ${\sigma }_i^2$. As in the previous case, the first step is the estimation of $(\mu ,{\mathit{\sigma }}^{\mathbf{2}})$ under the null hypothesis. In this case, we need to maximize (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ) under the restriction (22) ${\sigma }_1^2\ge {\sigma }_2^2\ge \cdots \ge {\sigma }_k^2>0.$(22)

Theorem 2

Suppose that for $1\le i\le k$, ${\overline{\sigma }}_i^2>0$. Then there exists a maximizer of (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ) under constraints (22(22) ${\sigma }_1^2\ge {\sigma }_2^2\ge \cdots \ge {\sigma }_k^2>0.$(22) ).

Proof.

First consider the situation for fixed ${\mathit{\sigma }}^{\mathbf{2}}$ with ${\sigma }_i^2>0$ for all i. Differentiating (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ) with respect to μ yields the equation $\sum _{i=1}^k\frac{n_i(\overline{y_i}-\mu )}{{\sigma }_i^2}$ This shows, that for this ${\mathit{\sigma }}^{\mathbf{2}}$, the (unique) maximizer of (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ) in μ is given by the following weighted sum of level-means, $\hat{\mu }\left({\mathit{\sigma }}^{\mathbf{2}}\right)=\frac{\sum _{i=1}^k{n}_i{\sigma }_i^{-2}{\overline{y}}_i}{\sum _{i=1}^k{n}_i{\sigma }_i^{-2}}$ Consequently, $\underset{i}{min}{\overline{y}}_i\le \hat{\mu }\left({\mathit{\sigma }}^{\mathbf{2}}\right)\le \underset{i}{max}{\overline{y}}_i$, bounding the set of possible maximizers of (18(18) $l(\mu ,{\mathit{\sigma }}^{\mathbf{2}})=\sum _{i=1}^k\left\{-\frac{n_i}{2}\mathrm{l}\mathrm{n}{\sigma }_i^2-\frac{1}{2{\sigma }_i^2}\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2\right\}+\mathrm{c}.$(18) ) in μ irrespective of the precise value of ${\mathit{\sigma }}^{\mathbf{2}}$.

Now, given any $\mu \in \mathbb{R}$, the corresponding optimal ${\mathit{\sigma }}^{\mathbf{2}}$ is the solution to the antitonic regression problem $\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{g}({\overline{\mathit{\sigma }}}_{{\mathit{H}}_{\mathbf{0}}}^{\mathbf{2}},\mathit{N})$ where ${\overline{\mathit{\sigma }}}^{\mathbf{2}}=({\overline{\sigma }}_1^2,\dots {\overline{\sigma }}_k^2{)}^{\prime }$, ${\overline{\sigma }}_i^2=\frac{\sum _{j=1}^{n_i}(y_{ij}-\mu {)}^2}{n_i}$, $\mathit{N}=(n_1,\dots ,n_k{)}^{\prime }$ (see [33, Example 1.5.5]). The vector to be projected has elements ${\overline{\sigma }}_i^2+(\mu -{\overline{y}}_i{)}^2$. This means, that if μ is restricted to $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$, the coordinates to be projected all belong to the interval $[\underset{i}{min}{\overline{\sigma }}_i^2,max{\overline{\sigma }}_i^2+(max{\overline{y}}_i-min{\overline{y}}_i{)}^2]$. So, if μ ranges over $[\underset{i}{min}{\overline{y}}_i,\underset{i}{max}{\overline{y}}_i]$, the optimal ${\mathit{\sigma }}^{\mathbf{2}}$ is also contained in a the closed bounded region $[\underset{i}{min}{\overline{\sigma }}_i^2,max{\overline{\sigma }}_i^2+(max{\overline{y}}_i-min{\overline{y}}_i{)}^2{]}^k$. By our assumption that all ${\overline{\sigma }}_i^2>0$, the MLE exists being a maximizer of a continuous function on a compact set in $\mathbb{R}\times {\mathbb{R}}^k$

If we consider this case as a special case of the case considered in [36], the solution is unique if Condition 2.2 holds. Given that the solution is not in a closed form, we use an iterative procedure to approximate the solution. As a starting value ${\hat{\mu }}^{\left(0\right)}$, a modified version of the Graybill-Deal estimator of the common mean when the variances are subject to order restrictions proposed in [25] appears to be a good choice: (23) ${\hat{\mu }}_{\left(I\right)}=\frac{\sum _{i=1}^k{w}_i{\hat{\tau }}_i{\overline{y}}_i}{\sum _{i=1}^k{w}_i{\hat{\tau }}_i}$(23) where ${\hat{\tau }}_i$ is the isotonic regression of $(\mathit{t},\mathit{N})$ where $\mathit{t}=(t_1,\dots t_k{)}^{\prime }$, $t_i=\frac{1}{s_i^2}$.

Under $H_1$ we use as estimates of $({\mathit{\mu }}_{{\mathit{H}}_{\mathbf{1}}},{\mathit{\sigma }}_{\mathit{i}\mathit{H}\mathbf{1}}^{\mathbf{2}})$, $({\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}^{\mathit{I}},{\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}^{\mathbf{2}\mathit{I}})$ found using the iterative procedure described in Section 2.3.

In contrast with the previous cases, it is not possible to further reduce the expression of the LRT because $\exp \left\{\frac{1}{2}\sum _{i=1}^k\sum _{j=1}^{n_i}\frac{(y_{ij}-{\hat{\mu }}_{H_0}{)}^2}{{\sigma }_{H_0}^2}\right\}$ does not reduce to a constant. The same holds under $H_1$. Therefore, the LRT in this case can be computed by substituting the solutions obtained via the iterative procedure under $H_0$ and $H_1$ in the generic expression given in (10(10) $\mathrm{\Lambda }=\frac{\underset{(\mathit{\mu }\in H_0;{\mathit{\sigma }}^{\mathbf{2}})}{max}L({\mathit{y}}_{\mathbf{1}},{\mathit{y}}_{\mathbf{2}},\dots ,{\mathit{y}}_{\mathit{k}};\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})}{\underset{(\mathit{\mu }\in H_1;{\mathit{\sigma }}^{\mathbf{2}})}{max}L({\mathit{y}}_{\mathbf{1}},{\mathit{y}}_{\mathbf{2}},\dots ,{\mathit{y}}_{\mathit{k}};\mathit{\mu },{\mathit{\sigma }}^{\mathbf{2}})}$(10) ): (24) ${\mathrm{\Lambda }}^I=\frac{L({\hat{\mu }}_{H_0},{\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}_{{\mathit{H}}_{\mathbf{0}}}^{\mathbf{2}\mathit{I}})}{L({\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}_{{\mathit{H}}_{\mathbf{1}}}^{\mathit{I}},{\stackrel{\mathbf{ˆ}}{\mathit{\sigma }}}_{{\mathit{H}}_{\mathbf{1}}}^{\mathbf{2}\mathit{I}})}$(24)

4. Null hypothesis distribution of the test statistics: bootstrap approach

In order to determine the significance of the various test statistics proposed in the previous sections, we need the null hypothesis distribution of the test statistics. The main distributional results concerning ${\overline{\chi }}_k^2$ and ${\overline{E}}_k^2$, the test statistics derived in the known variance ratio case, are contained in [2, Theorems 3.1–3.2]. However, problems related to the value of k can arise in the analytical derivation of the p-values. Numerical approximation can be necessary, especially if k>4 and if the variation in the range of the weights is not ‘moderate’ [32,39].

Furthermore, in the case of completely unknown variances, the null distribution depends on the unknown variances. When analytical derivation of the null distribution is particularly complex or not possible, bootstrap methodology is a good option. Therefore, we propose both a parametric and a non-parametric bootstrap approach that can be easily employed for finding approximate p-values taking into account the different assumptions on the variances. For overcoming the complex derivation when the variances are unknown, bootstrap procedures have been proposed in the literature [4,24].

In particular, in [24], an interesting review of the methods used to approximate the null distribution of the test statistic under $H_0$ and the restrictive normality assumption (with which we will not deal in this paper) is reported. Moreover, the authors propose both a parametric and non-parametric bootstrap approach for the LRT null distribution for one-sided hypothesis testing for means in a multivariate setting [24]. Also in [4], a bootstrap approach to test the homogeneity of order restricted mean vectors when the covariance matrices are unknown is used. In line with those previous approaches, here we propose two general bootstrap procedures, parametric and non-parametric, that can be used for testing the null hypothesis taking into account the various assumptions on the variances.

4.0.1. Parametric bootstrap

4.0.2. Non-parametric bootstrap

The non-parametric version of the bootstrap releases the normality assumption of the bootstrap samples. However, a relatively large sample size is required for the following approach.

5. Application

One of the most common ways for investigating strength and ductility of metallic materials is by performing a tensile test. Loosely speaking, a tensile test is an experiment in which force is applied to the test sample causing deformation of the material, temporarily (elastic behavior), permanently (plastic behavior) and eventually its fracture [9]. Data used in this paper are image data of the microstructure of the material subjected to a plastic strain (deformation) of 0.139 obtained performing a uniaxial tensile test in which force is applied to the test sample with respect to just one specific axis (Figure 1).

Figure 1. Tensile testing machine.

At a microstructure level, the deformation of the material corresponds to displacements in the lattice structure and in the possible appearance of GNDs. The material used in this paper is an annealed AISI420 stainless steel with $M_{23}{C}_6$ carbides and aim is investigating the carbide effect on the GNDs formation. KAM is used as a proxy of the GNDs. In Figure 2, the KAM is represented by filaments inside and among grains delimited by solid lines called grain boundaries; carbides are the black dots.

Figure 2. Microstructure image showing the KAM at strain level 13.9% (overlapped grid of $25\times 25$).

Modeling the relationship between KAM and carbides and more generally understanding its inhomogeneous distribution over the microstructure is now the main aim and it can be considered a starting point for finding a stochastic model for predicting mechanical properties from 2D microstructure images. We apply estimation procedures and perform tests under order restrictions, solving the three different univariate isotonic regressions problems (2.1, 2.2, 2.3) according to the assumption on the variances.

The first step for obtaining the data in the most suitable form for the analysis is ‘overlaying’ a grid over the image. In Figure 2, a $25\times 25$ grid is added to the image. With $y_{ij}$, we denote the mean KAM value (expressed in misorientation degree $[0,5]$) of the jth square of the grid of the image taken in which i carbides are observed. The explanatory variable $\mathrm{X}$ in all three isotonic regressions is the number of carbides observed in the grid squares. A plot of the data is shown in Figure 3.

Figure 3. Plot of KAM and number of carbides for the 625 squared areas of Figure 2.

We wish to test the null hypothesis that the expected KAM is the same in all the squares of the gird, regardless the numbers of carbides observed in the grid. The alternative hypothesis (26) ${\mu }_0\le {\mu }_1\le {\mu }_2\le {\mu }_3.$(26) represents the idea that KAM tends to be higher in areas where more carbides are observed. Moreover, in the ordered variances case, we assume that (27) ${\sigma }_0^2\ge {\sigma }_1^2\ge {\sigma }_2^2\ge {\sigma }_3^2>0.$(27) This is in accordance with what we see in Figure 3. In fact, the idea behind this assumption is that in areas in which less carbides are observed GNDs have more freedom to move, resulting in increments in dispersion.

The results of the estimation in the three different scenarios faced in Sections 2.1–2.3 are summarized in Table 1. The values of $\overline{y}$, ${\overline{\sigma }}^2$ and $s^2$ represent the empirical group mean, group variance and group sample variance respectively. Comparing the solutions of the three different isotonic regression, (1) known variance ratio, (2) unknown variance and (3) variances under order restrictions, just very slightly differences can be noticed.

Table 1. Values of estimated means and variances of the KAM conditioned on the number of carbides visible in a square of a grid $25\times 25$ according to different order restrictions assumptions (13.9% Strain).

	0	1	2	3
$\overline{y}$	0.815	0.833	0.870	0.854
${\overline{\sigma }}^2$	0.035	0.024	0.017	0.022
$s^2$	0.035	0.024	0.017	0.023
${\hat{\mu }}_{\left(1\right)}^I$	0.815	0.833	0.867	0.867
${\hat{\mu }}_{\left(2\right)}^I$	0.815	0.833	0.867	0.867
${\hat{\sigma }}_{\left(2\right)}$	0.035	0.024	0.017	0.022
${\hat{\mu }}_{3)}^I$	0.815	0.833	0.866	0.866
${\hat{\sigma }}_{\left(3\right)}^I$	0.035	0.024	0.018	0.018
n	340	211	54	18

In Table 2, the results of testing the null hypothesis are shown. For computing ${\chi }^{2\left(0\right)}$ and ${\overline{E}}^{2\left(0\right)}$, the variance ratio is supposed to be known. In the specific case, we assume that the KAM total variance for the whole image is the real-known variance and that $c_i=\frac{{\overline{\sigma }}_i^2}{{\sigma }^2}$. For computing both the parametric and non-parametric p-values, the two different bootstrap approaches described in Section 4 have been used and M, the number of replications, is taken equal to 20, 000. Independent of the knowledge or assumptions on the variances, the conclusion is the same and leads to the rejection of the null hypothesis.

Table 2. Estimated values for the four different LRT with the corresponding parametric and non-parametric p-values.

		Test	p-value	p-value
	${\hat{\mu }}_{H_0}$	statistic	(parametric)	(non-parametric)
${\overline{\chi }}^{2\left(0\right)}$	0.827	5.760	0.0323	0.0310
${\overline{E}}^{2\left(0\right)}$	0.831	0.0121	0.0112	0.0085
$-2\log \tilde{\mathrm{\Lambda }}$	0.831	7.330	0.0178	0.0222
${\mathrm{\Lambda }}^{I\left(0\right)}$	0.831	7.105	0.0212	0.0251
	0	1	2	3
${\hat{\sigma }}_{H_0}^2$	0.035	0.024	0.018	0.022
${\hat{\sigma }}_{H_0}^{2I}$	0.035	0.024	0.019	0.019

6. Conclusions

This paper presents three different models involving order restrictions and within these models the ML estimators and LRTs for the homogeneity of the means against monotonicity are introduced and studied. Prior knowledge given by physical relations or intuition is not often exploited in statistical studies about materials and this can lead to less efficient methods that produce less accurate results. After having described the estimation procedures and highlighted how prior knowledge of the variances influence these, we propose the LRT as test statistic for testing homogeneity of means. In the unknown variances case and the ordered variances case, heteroskedasticity plays a crucial role also under the null hypothesis, leading to different estimates of the common mean under $H_0$. Results on existence and uniqueness of the maximum likelihood estimates in these last two cases are derived. Furthermore, two different bootstrap approaches are proposed for approximating the null distribution of the test statistic under the different assumptions on the variances. The proposed tests are applied to a real data example from Materials Science, showing evidence that the so-called KAM tends to be higher in regions of the microstructure where more carbides are observed. In fact, incorporating reasonable intuition about the order of means and the variances order in this context helps understanding the evolution of complicated structure of dislocations in metals and its effect on the hardening behavior of the material during deformation.

Acknowledgments

This research was carried out under project number S41.5.14547b in the framework of the Partnership Program of the Materials innovation institute M2i (www.m2i.nl) and the Technology Foundation TTW (www.stw.nl), which is part of the Netherlands Organization for Scientific Research (www.nwo.nl).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was carried out under project number S41.5.14547b in the framework of the Partnership Program of the Materials innovation institute M2i (www.m2i.nl) and the Technology Foundation TTW (www.stw.nl), which is part of the Netherlands Organization for Scientific Research (www.nwo.nl).

Notes

1 Dislocations are usually classified into redundant and non-redundant dislocations, respectively, called Statistically Stored Dislocations (SSDs) and Geometrically Necessary Dislocations (GNDs). GNDs are dislocations with a cumulative effect and they allow the accommodation of lattice curvature due to non-homogeneous deformation. They control the work hardening individually by acting as obstacles to slip and collectively by creating a long-range back stress.

Appendix 1

A.1. ALGORITHM 2.1

• Let ${\mathit{\mu }}^{\left(0\right)}=\overline{\mathit{y}}=({\overline{y}}_1,\dots {\overline{y}}_k{)}^{\prime }$ and $\mathit{w}=(w_1,\dots w_k{)}^{\prime }$, $w_i=\frac{n_i}{c_i}$

• Is ${\overline{y}}_1\le {\overline{y}}_2\le \cdots \le {\overline{y}}_k$?

  1. ${\mathit{\mu }}^{\ast }=\overline{\mathit{y}}$ is the solution

  2. ${\overline{y}}_i>{\overline{y}}_{i+1}$

     Replace ${\overline{y}}_i$ and ${\overline{y}}_{i+1}$ by (A1) $m_{i,i+1}=\frac{w_i{\overline{y}}_i+w_{i+1}{\overline{y}}_{i+1}}{w_i+w_{i+1}}$(A1) Repeat until QUESTION 1 is satisfied.

A.2. ALGORITHM 2.2

A.2.0.1 Two steps Iterative procedure

• Let ${\mathit{\mu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=\overline{\mathit{y}}=({\overline{y}}_1,\dots {\overline{y}}_k{)}^{\prime }$, ${\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathbf{0}\mathbf{\right)}}={\overline{\mathit{\sigma }}}^{\mathbf{2}}=({\overline{\sigma }}_1^2,\dots {\overline{\sigma }}_k^2{)}^{\prime }$, ${\overline{\sigma }}_i^2=\frac{\sum _j(y_{ij}-{\overline{y}}_i{)}^2}{n_i}$ and ${\mathit{w}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=(w_1^{\left(0\right)},\dots w_k^{\left(0\right)}{)}^{\prime }$, $w_i^{\left(0\right)}=\frac{n_i}{{\sigma }^{2\left(0\right)}}$

• Is ${\overline{y}}_1\le {\overline{y}}_2\le \cdots \le {\overline{y}}_k$?

  1. ${\mathit{\mu }}^{\mathbf{\ast }}={\mathit{\mu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=\overline{\mathit{y}}$ and ${\mathit{\sigma }}^{\mathbf{2}}={\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathbf{0}\mathbf{\right)}}={\overline{\mathit{\sigma }}}^{\mathbf{2}}$ are the solutions

  2. Use Step 1.2 Algorithm 2.1 to compute ${\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ with weights ${\mathit{w}}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}}$

  3. Compute ${\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathit{l}\mathbf{\right)}}={\mathit{s}}^{\mathbf{2}\mathbf{\left(}\mathit{l}\mathbf{\right)}}$, $s_i^{2\left(l\right)}=\frac{\sum _j(y_{ij}-{\mu }_i^{\left(l\right)}{)}^2}{n_i}$ and ${\mathit{w}}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}=(w_1^{\left(l\right)},\dots w_k^{\left(l\right)}{)}^{\prime }$, $w_i^{\left(l\right)}=\frac{n_i}{{\sigma }^{2\left(l\right)}}$

  4. Go back to QUESTION 1 using ${\mathit{w}}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$.

     Repeat until $\underset{1\le i\le k}{max}|{\mu }_i^{\ast (l-1)}-{\mu }_i^{\ast \left(l\right)}|\le {10}^{-m}$

A.2.0.2 Alternating Iterative Method

• Let ${\mathit{\nu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=(1/{\overline{\sigma }}_1^2,\dots 1/{\overline{\sigma }}_k^2{)}^{\prime }$, ${\overline{\sigma }}_i^2=\frac{\sum _j(y_{ij}-{\overline{y}}_i{)}^2}{n_i}$

• ${\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ the isotonic regression on $D_a$ using weights ${\mathit{w}}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}}=(w_1^{(l-1)},\dots w_k^{(l-1)}{)}^{\prime }$, $w_i^{(l-1)}=n_i{\nu }^{(l-1)}$;

• ${\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ maximizing the profile likelihood $L(y;{\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}};\mathit{\nu })$ on $V_0$, $V_0=\{\nu \in {\mathbb{R}}^k:0\le 1/\underset{i}{max}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right))\le {\nu }_1\le \cdots \le {\nu }_k\le 1/\underset{i}{min}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)\left)\right\}$.

  Repeat (1)-(2) until $|L(y;{\mathit{\mu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}},{\mathit{\nu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}})-L(y;{\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}},{\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}})|\le {10}^{-m}$

A.3. ALGORITHM 2.3

A.3.0.3 Two steps Iterative procedure

• Let ${\mathit{\mu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=\overline{\mathit{y}}=({\overline{y}}_1,\dots {\overline{y}}_k{)}^{\prime }$, ${\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathbf{0}\mathbf{\right)}}={\overline{\mathit{\sigma }}}^{\mathbf{2}}=({\overline{\sigma }}_1^2,\dots {\overline{\sigma }}_k^2{)}^{\prime }$, ${\overline{\sigma }}_i^2=\frac{\sum _j(y_{ij}-{\overline{y}}_i{)}^2}{n_i}$ and and ${\mathit{w}}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=(w_1^{\left(0\right)},\dots w_k^{\left(0\right)}{)}^{\prime }$, $w_i^{\left(0\right)}=\frac{n_i}{{\sigma }^{2\left(0\right)}}$.

• Is ${\overline{y}}_1\le {\overline{y}}_2\le \cdots \le {\overline{y}}_k$?

  1. ${\mathit{\mu }}^{\mathbf{\ast }}={\mathit{\mu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=\overline{\mathit{y}}$ go to QUESTION 2.

  2. Use Step 1.2 Algorithm 2.1 to compute ${\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ with weights ${\mathit{w}}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}}$

  3. Compute ${\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathit{l}\mathbf{\right)}}={\mathit{s}}^{\mathbf{2}\mathbf{\left(}\mathit{l}\mathbf{\right)}}$, $s_i^{2\left(l\right)}=\frac{\sum _j(y_{ij}-{\mu }_i^{\left(l\right)}{)}^2}{n_i}$ and ${\mathit{w}}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}=(w_1^{\left(l\right)},\dots w_k^{\left(l\right)}{)}^{\prime }$, $w_i^{\left(l\right)}=\frac{n_i}{{\sigma }^{2\left(l\right)}}$

• Is ${\sigma }_1^{2\left(l\right)}\ge {\sigma }_2^{2\left(l\right)}\ge \cdots \ge {\sigma }_k^{2\left(l\right)}$?

  1. ${\mathit{\mu }}^{\mathbf{\ast }}={\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ and ${\mathit{\sigma }}^{\mathbf{2}\mathbf{\ast }}={\mathit{\sigma }}^{\mathbf{2}\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ are the solutions.

  2. ${\sigma }_i^{2\left(l\right)}<{\sigma }_j^{2\left(l\right)}$

  3. Replace ${\sigma }_i^{2\left(l\right)}$ and ${\sigma }_{i+1}^{2\left(l\right)}$ by (A2) ${\sigma }_i^{2(l+1)}={\sigma }_{i+1}^{2(l+1)}={\overline{s}}_{i,i+1}=\frac{n_i{\sigma }_i^{2\left(l\right)}+n_{i+1}{\sigma }_{i+1}^{2\left(l\right)}}{n_i+n_{i+1}}$(A2) Repeat until QUESTION 2 is satisfied.

  4. Go back to QUESTION 1 using ${\mathit{w}}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$.

     Repeat until $\underset{1\le i\le k}{max}|{\mu }_i^{\ast (l-1)}-{\mu }_i^{\ast \left(l\right)}|\le {10}^{-m}\mathrm{a}\mathrm{n}\mathrm{d}\underset{1\le i\le k}{max}|{\sigma }_i^{2\ast (l-1)}-{\sigma }_i^{2\ast \left(l\right)}|\le {10}^{-m}$

A.3.0.4 Alternating Iterative Method

• Let ${\mathit{\nu }}^{\mathbf{\left(}\mathbf{0}\mathbf{\right)}}=(1/{\overline{\sigma }}_1^2,\dots 1/{\overline{\sigma }}_k^2{)}^{\prime }$, ${\overline{\sigma }}_i^2=\frac{\sum _j(y_{ij}-{\overline{y}}_i{)}^2}{n_i}$

• ${\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ use Step 1 AlM Algorithm 2.2;

• ${\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}}$ the isotonic regression on $V_0$, $V_0=\{\nu \in {\mathbb{R}}^k:0\le 1/\underset{i}{max}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right))\le {\nu }_1\le \cdots \le {\nu }_k\le 1/\underset{i}{min}(\underset{min\left(\overline{\mathit{y}}\right)\le \theta \le max\left(\overline{\mathit{y}}\right)}{min}{s}_i^2\left(\theta \right)\left)\right\}$ with weights $\mathit{N}=(n_1,\dots ,n_k{)}^{\prime }$ Repeat (1)–(2) until $|L(y;{\mathit{\mu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}},{\mathit{\nu }}^{\mathbf{(}\mathit{l}\mathbf{-}\mathbf{1}\mathbf{)}})-L(y;{\mathit{\mu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}},{\mathit{\nu }}^{\mathbf{\left(}\mathit{l}\mathbf{\right)}})|\le {10}^{-m}$

Arxiv Preprint: arXiv:2002.00993