An Adomian decomposition method with some orthogonal polynomials to solve nonhomogeneous fractional differential equations (FDEs)

1.   Introduction

In practice, many engineering or medical studies are concerned with comparing different group means against each other or grand mean. Several parametric and nonparametric methods are used for this purpose ^[1]. The analysis of variance (ANOVA) technique based on the F-statistic is a well-known one-way fixed effect method to differentiate group means ^[2]. In addition, there are various statistical procedures utilized for comparing independent group means, such as, the Welch test ^[3], the James-second-order test ^[4], Brown-Forsythe test ^[5], and Alexander-Govern test ^[6]. The structure of the preceding tests is compatible with assessing the pairwise significance of treatment means. An extension of ANOVA in the presence of covariate is named analysis of covariance (ANCOVA), which is used to examine whether there is a statistically significant difference between the means of three or more independent groups after taking into account one or more covariates ^[7,8,9]. However, for examining the difference in treatment means from their grand mean, we use the analysis of means (ANOM) test originated and formally proposed by Ott ^[10] (reproduced by Ott ^[11]). The ANOM test is applied for analysis in several fields, such as environmental studies ^[12], medical science ^[13,14], nanomaterials ^[15], tourism ^[16], and healthcare studies ^[17]. The extension of the ANOM test under mixed effect designs and balanced incomplete block designs was proposed by Schilling ^[18]. The ANOM test is a graphical method that is not only useful for comparing group means but also beneficial for comparing rates or proportions ^[19].

Initially, the ANOM test was designed for the equality of means; Wludyka and Nelson ^[20] proposed the ANOM mechanism for the equality of variances, which is known as the analysis of means for variances (ANOMV). Bernard and Wludyka ^[21] and Wludyka and Sa ^[22] suggested the robustness of ANOMV with the combination of the Fligner and Killeen test and the Levene test. An extension of the ANOM test under a heteroscedastic model, known as heteroscedastic analysis of means (HANOM), was proposed by Nelson and Dudewicz ^[23] and Dudewicz and Nelson ^[24]. A nonparametric version of the ANOM test was introduced by Bakir ^[25], and a comparison between ANOM and ANOVA tests using parametric bootstrap was conducted by Chang et al. ^[26] The exact control limits for the balanced design with equal sample sizes were presented by Nelson ^[27], Nelson ^[28], while for the unbalanced design with unequal sample sizes were given by Soong and Hsu ^[29]. Further, the tables for the ANOM test with equal sample sizes were reported in studies ^{[30,31,32,33]} and for unequal sample sizes in studies ^[34,35].

Recently, Mendeş and Yiğit ^[36] established a comparative study between ANOVA-F and ANOM tests under the violation of assumptions (e.g., normality, homogeneity of variances) in terms of type Ⅰ error rate and power of the test. Guirguis and Tobias ^[37] produced the distributional properties of the ANOM test using Fortran syntax, and Pallmann and Hothorn ^[38] presented the applications of the ANOM test by using the R language. ^[38] introduced the generalized approach for ANOM utilizing the concept of multiple contrasts tests (MCTs), specific comparisons to the grand mean, and further generalizations for MCTs by using a linear model with a covariate. The previous ANOM versions were considered for the fixed-effect model, while Jayalath and Ng ^[39] examined the ANOM test for the random effect model, and Jayalath and Ng ^[40] proposed the ANOM test for hierarchically nested and split-plot designs. A brief literature review on the ANOM test can be found in ^[41]. An individual measurement control chart based on ANOM control limits was suggested by Chakravarthi and Rao ^[42]. The effect of measurement errors on the performance of the ANOM test was studied by Chakraborty and Khurshid ^[43], and a Bayesian graphical approach for the location parameter of the process was discussed by Apley ^[44]. The bootstraps confidence interval of the ANOM and ANOVA were derived by Lopez-Mejia and Roldan-Valadez ^[45].

Generally, many experiments contain a study variable $\left(Y\right)$ that is observed with another linearly associated variable $\left(X\right)$. The variable $X$ is known as a covariate or concomitant variable, which is an uncontrollable predictor and is found along with the study variable ^{[46,47,48,49,50,51,52,53]}. These types of variables are common in many fields, such as: in the monofilament fibre or glue industry; the strength (study variable) produced by different machines is affected by the thickness (concomitant variable) of the fiber in the cutting machine; the amount of metal removed (study variable) is associated with the hardness of the specimen (concomitant variable), in Electrical Discharge Machining (EDM) problem, the effect of pulse current on the surface roughness (study variable) is affected by the machining time (concomitant variable), in medical science; effect of Viagra dosage on participants libido (study variable) is affected by the partners' libido (concomitant variable) and in marine studies; growth (weight) of oyster (study variable) is dependent on the initial weight of oyster (concomitant variable) ^{[54,55,56,57]}.

From the above-stated literature, it can be seen that the traditional ANOM test does not consider the concomitant variable that may disturb the mean square error and, consequently, may conduct false judgments about the potential differences among different treatments. In this study, we intend to propose a new testing mechanism named the analysis of means with covariate (ANOMC). The new technique is developed under the following scenarios:

1) Measure the study variable $\left(Y\right)$ and a covariate $\left(X\right)$ among several groups (or treatments).

2) Assume a linear relationship between $Y$ and $X$ for each group.

3) Compare treatment adjusted means against their grand mean conditional on the value of $X$.

4) Identify which treatment's adjusted mean is exactly significant.

The newly proposed methodology will give an indication of the significant mean using adjusted mean effects.

The rest of the article is organized as follows. In Section 2, we describe the brief methodology of ANOM and ANOMC tests. The design parameters of the study are reported in Section 3. Section 4 evaluates the performance of the proposed and competing methods. Section 5 presents illustrative examples of mechanical/industrial engineering and medical phenomena. Finally, Section 6 provides a summary, conclusions, and recommendations for the study.

2.   Description of existing and proposed methods

In this section, firstly, we will outline the methodology of the traditional ANOM test about testing of means. Later, we will describe the newly proposed method named by the ANOMC test for the testing of adjusted means in the presence of a covariate.

2.1.   The analysis of means (ANOM)

For the completely randomized design (CRD) with a single-factor model having $t$ treatments, each with ${n}_{i}$ observations, and the total number of observations is $N = \sum _{i = 1}^{t}{n}_{i}$. The fixed-effects model can be represented as follows:

where ${Y}_{ij}$ is the ${{\mathrm{j}}}^{th}$ observation of response variable for the ${i}^{th}$ treatment level, ${\tau }_{i} = \mu -{\mu }_{i}$ is the fixed effect of the ${i}^{th}$ treatment level from overall mean $\mu$. ${ϵ}_{ij}$ is the ${{\mathrm{j}}}^{th}$ random error of the ${i}^{th}$ treatment level and assumed to be normally distributed with zero mean and constant variance ${\sigma }^{2}$. The variance ${\sigma }^{2}$ is assumed to be constant for all treatment levels, which implies that the observations ${Y}_{ij}~N(\mu +{\tau }_{i}, {\sigma }^{2})$ and the observations are mutually independent. The model given in Eq (2.1) is a statistical linear model, i.e., the response variable ${Y}_{ij}$ is a linear function of the model parameters. The layout of the ANOM data set is presented in Table 1, where the structure of the ANOM test under the same assumptions as the model in Eq (2.1) is used to test the following hypotheses:

Null hypothesis, ${H}_{0}:{\mu }_{1} = {\mu }_{2} = \dots = {\mu }_{t}$;

Alternative hypothesis, ${H}_{1}:$ at least one group mean differs from the grand mean.

Under the balanced design (equal sample sizes among all treatments $\left({n}_{i} = n\right)$), the lower decision line $\left(LDL\right)$ and upper decision line $\left(UDL\right)$ for ANOM test are defined as below:

where ${\overline{Y}}_{..} = \sum _{i = 1}^{t}{\overline{Y}}_{i.}/t$ is the grand mean, $MSE = \sum _{i = 1}^{t}{S}_{{Y}_{i.}}^{2}/t$ is the mean square error, $\alpha$ is the pre-specified type Ⅰ error rate, $n$ is the sample size, ${n}_{i}$ is the sample size of ${i}^{th}$ treatment, $N = \sum _{i = 1}^{t}{n}_{i}$ is the total number of observations, $t$ is the number of treatments, and $h(\alpha, t, N-t)$ is the critical value reported in Table B.1 ^[58].

However, under the unbalanced design (unequal sample sizes), the lower decision line $\left(LDL\right)$ and upper decision line $\left(UDL\right)$ for the ANOM test are expressed as below:

where ${\overline{Y}}_{..} = {n}_{1}{\overline{Y}}_{1.}+{n}_{2}{\overline{Y}}_{2.}+\dots +{n}_{t}{\overline{Y}}_{t.}/N$ is the weighted overall mean, $MSE = \left({n}_{1}-1\right){S}_{{Y}_{1.}}^{2}+\left({n}_{2}-1\right){S}_{{Y}_{2.}}^{2}+\dots +({n}_{t}-1){S}_{{Y}_{t.}}^{2}/N$ is the pooled mean square error, and $m(\alpha, t, N-t)$ is the critical value reported in Table B.3 ^[58], and all other notations are the same as discussed above. The ${i}^{th}$ treatment mean is declared significantly different from the grand mean if ${\overline{Y}}_{i.}$ falls outside the $LDL$ and $UDL$.

It is to be mentioned that the ANOM test assumes normality and homogeneity of the variances, which are briefly discussed and compared with the ANOVA-F test by Mendeş and Yiğit ^[36]. It is to be noted that Tables B.1 and B.3 given in Nelson et al. ^[58] have critical values with limited choices of the degree of freedom. So, by adopting Nelson ^[28] mechanism, we have derived the critical values $h\left(\alpha, t, N-t\right)$ and $m(\alpha, t, N-t)$ (given in Table 2) with the parameter choices considered in this study.

2.2.   The analysis of means with covariate (ANOMC)

Assume a single-factor model with a linearly related covariate having $T$ treatments with $t$ levels, each with ${n}_{i}$ observations, and the total number of observations is $N = \sum _{i = 1}^{t}{n}_{i}$. The fixed-effects model can be represented as follows:

where ${Y}_{ij}$ is the ${{\mathrm{j}}}^{th}$ observation of response variable for the ${i}^{th}$ treatment level, ${X}_{ij}$ is the ${{\mathrm{j}}}^{th}$ observation of covariate for the ${i}^{th}$ treatment level corresponding to ${Y}_{ij}$. Further, $\mu$ is the overall mean, ${\tau }_{i}$ is the effect of ${i}^{th}$ treatment level, $B$ is the slope indicating the relationship between ${Y}_{ij}$ and ${X}_{ij}$, ${\overline{X}}_{..}$ is the mean of ${X}_{ij}$ observations and ${ϵ}_{ij}$ is the ${{\mathrm{j}}}^{th}$ random error of the ${i}^{th}$ treatment level and assumed to be normally distributed with zero mean and constant variance ${\sigma }^{2}$. It is noted that in the model (2.6), we assumed that the slope $B\ne 0$, and the relationship between ${Y}_{ij}$ and ${X}_{ij}$ is linear, the regression coefficients for each treatment are identical, the concomitant variable ${X}_{ij}$ is not affected by treatment, and the treatment effects sum to zero (i.e., ${\sum }_{i = 1}^{t}{\tau }_{i} = 0)$. The model given in Eq (2.6) is also a statistical linear model, i.e., the response variable ${Y}_{ij}$ is a linear function of the model parameters.

The layout of the ANOMC dataset is reported in Table 3. The ANOMC test under the same assumptions as the model in Eq (2.6) is used to test the following hypotheses:

Null hypothesis ${, H}_{0}:{\mu }_{1} = {\mu }_{2} = \dots = {\mu }_{t}$;

Alternative hypothesis, ${H}_{1}:$ at least one group mean differs from the grand mean.

In the ANOMC test, the adjusted $\left(Adj\right)$ means are calculated by the regression mean estimator $\left({\overline{M}}_{i.}\right),$ which is an unbiased estimator (i.e., ${\overline{M}}_{..} = \sum _{i = 1}^{t}{\overline{M}}_{i.}/t = {\overline{Y}}_{..}$) having the minimum standard deviation $\left({S}_{{\overline{M}}_{i.}}\right)$. For more details about regression estimators, see ^{[59,60,61,62,63]}.

Under the balanced design (equal sample sizes), the lower decision line $\left(LDL\right)$ and upper decision line $\left(UDL\right)$ for the ANOMC test are defined below:

where ${h}^{{\mathrm{*}}}(\alpha, t, n)$ is the critical value reported in Table 2 (for more details, see Section 3.1). The ${b}_{i.}$ is the ${i}^{th}$ slope, ${\overline{M}}_{i.}$ is the ${i}^{th}$ regression mean estimator, ${r}_{i.}$ is the ${i}^{th}$ sample correlation, ${k}_{i.}$ is the ${i}^{th}$ unbiasing constant, ${S}_{{\overline{M}}_{i.}}$ is the sample standard deviation of ${i}^{th}$ regression mean estimator, ${\overline{M}}_{..} = \sum _{i = 1}^{t}{\overline{M}}_{i.}/t$ is the grand regression mean and ${MSE}_{{\overline{M}}_{..}} = \sum _{i = 1}^{t}{S}_{{\overline{M}}_{i.}}/t$ is the overall mean square error.

Under the unbalanced design (unequal sample sizes), the lower decision line $\left(LDL\right)$ and upper decision line $\left(UDL\right)$ for the ANOMC test are calculated by the following expressions:

where ${m}^{{\mathrm{*}}}(\alpha, t)$ is the critical value reported in Table 2 (for more details, see Section 3.1). However, ${\overline{M}}_{..}$ and ${MSE}_{{\overline{M}}_{..}}$ are defined as follows:

It is to be noted that under the unbalanced design (unequal sample sizes), the values of the lower decision line $\left(LDL\right)$ and upper decision line $\left(UDL\right)$ vary with sample size. The ${i}^{th}$ treatment adjusted mean declared significantly different from the adjusted grand mean if ${\overline{M}}_{i.}$ falls outside the $LDL$ and $UDL$.

As mentioned earlier, the ANOM test worked under some assumptions such as normality, homogeneity, and linear relationship see ^[64,65,66]. Similarly, the ANOMC test also works under some assumptions, including (ⅰ) normal distribution of the study variable for each value of covariate variable within each treatment group, (ⅱ) variances of the conditional study variable are the same for each treatment group (Homogeneity), (ⅲ) linear relation between the study variable and covariate and (ⅳ) the regression coefficients for each treatment are identical (homogeneity of regression slopes).

3.   Design of the study

This section provides the design structure of this study, which is further considered to execute the simulation study. In this procedure, normal random numbers $\left({Z}_{j}; j = {\mathrm{1, 2}}, 3, \dots, n\right)$ having parameters $\left(\mu = 0;\sigma = 1\right)$ are generated by using the Box and Muller transformation ^[67,68]. The non-normal numbers are generated by the Flieshman mechanism ^[69] with four specified moments. Flieshman method defines a random number by using the polynomial transformation equation, which is given as follows:

where $a, b, c,$ and $d$ are the coefficients of transformation (cf. Table 4), and $V$ is the resulting variable having zero mean, unit variance with specified skewness $\left(s\right)$ and kurtosis $\left(kr\right)$ values. The Flieshman's transformation coefficients for the specified pair of skewness and kurtosis $(s, kr)$ are reported in Table 2. The skewness and kurtosis pair $\left(s, kr\right)$ is used to describe different distributions having zero mean and unit variance, such as $(0, \; \; \; \; 0)$ forms normal distribution, $(0, \; \; \; \; 3)$ forms heavy-tailed double exponential distribution, $(2, \; 6)$ forms extremely positive skewed exponential distribution, and $(0, \; \; \; \; 25)$ forms very heavy-tailed approximately Cauchy distribution. To get the desired slope and homogeneous values of a concomitant variable, the following model is used:

where ${Y}_{tj}$ is the ${j}^{th}$ response observation of ${t}^{th}$ treatment, ${\rho }_{t}$ is the correlation between $Y$ and $X$ for ${t}^{th}$ treatment, ${X}_{tj}$ is the ${j}^{th}$ observation of concomitant variable associated with ${t}^{th}$ treatment and ${E}_{tj}$ is the ${j}^{th}$ observation of error term associated with ${t}^{th}$ treatment. It is noted that both ${X}_{tj}$ and ${E}_{tj}$ are obtained by using the algorithm of Eq (3.1).

Under the procedural description given above, we have assessed the performance of ANOMC and ANOM tests under several aspects, including the following:

1) normality; different choices of $D = (s, kr)$,

2) correlation; different choices of correlation between $Y$ and $X\left(\rho \right)$,

3) homogeneity of variances; several cases of variances $\left(v\right)$,

4) hypotheses; null case $\left({\delta }_{1}\right)$ and non-null cases $\left({\delta }_{2-10}\right)$,

5) number of treatments $\left(t\right)$,

6) sample size $\left(n\right)$.

Further, the choices of the aforementioned design parameters $\left(\left(s, kr\right), \rho, v, \delta, t{\mathrm{a}}{\mathrm{n}}{\mathrm{d}}n\right)$ are presented in Table 5. In Table 5, the symbol ":" is used to differentiate the value of each treatment. For example, ${\rho }_{2} = 0.5:0.5:0.8:0.8$ means that the correlation between $Y$ and $X$ in the first two treatments is 0.5, and in the last two treatments, it is set at 0.8. It is to be noted that only one covariate is used in this study, and its distribution is assumed to be standard normal throughout the study. The Monte Carlo simulation study (motivated by ^[70]) is carried out by using the R software version (4.0.3).

As mentioned above, the decision lines (i.e., LDL and UDL) of the ANOMC method depend on the critical values. The procedure to find the critical values for the ANOMC method is illustrated in the following steps:

1) On the fixed correlation $\left({\rho }_{1}\right)$, variance ratio $\left({v}_{1}\right)$ and pair of skewness and kurtosis $\left({D}_{1} = \left({\mathrm{0, 0}}\right)\right)$, choose any case of the number of treatments $\left(t\right)$ and sample size $\left(n\right)$, under the null hypotheses.

2) Generate random numbers based on the information assumed in the previous step using the Flieshman method.

3) Calculate the statistics ${\overline{M}}_{i.}$, ${\overline{M}}_{..}$ and ${MSE}_{{\overline{M}}_{..}}.$

4) Use an arbitrary value as a critical value (i.e., ${h}^{{\mathrm{*}}}(\alpha, t, n)$ for balanced design and ${m}^{{\mathrm{*}}}\left(\alpha, t\right)$ for unbalanced design), and obtain a lower decision line $\left(LDL\right)$ and an upper decision line $\left(UDL\right)$ for the ANOMC test.

5) Plot the ${\overline{M}}_{i.}$ against the decision lines. Further, calculate an indicator variable $I$ such that it can have an observation equal to one, if any ${\overline{M}}_{i.}$ falls outside of the decision lines; otherwise, assumed to equal zero.

6) Repeat steps 1–5, a large number of runs to obtain specified $\alpha$.

If specified $\alpha$ does not achieve, then adjust the previous arbitrary critical value and repeat steps 1–6 until specified $\alpha$ is obtained. The obtained critical values for the ANOMC test are reported in Table 2 with respect to $\alpha = 0.01$ and $\alpha = 0.05$.

4.   Performance analysis

The performance of the two methods is investigated in terms of percentage type Ⅰ error $\left(\alpha \right)$ and the percentage power of the test $\left(1-\beta \right)$ ^[71]. The type Ⅰ error is the degree of the incorrect rejection of a true null hypothesis $\left({H}_{0}:{\mu }_{1} = {\mu }_{2} = \dots { = \mu }_{t}\right)$ which is mathematically defined as:

However, the power of the test is the degree of correct rejection of the false null hypothesis $\left({H}_{1}:\; {\mathrm{a}}{\mathrm{t}}\; {\mathrm{l}}{\mathrm{e}}{\mathrm{a}}{\mathrm{s}}{\mathrm{t}}\; {\mathrm{o}}{\mathrm{n}}{\mathrm{e}}\; {\mathrm{o}}{\mathrm{f}}\; {\mathrm{t}}{\mathrm{h}}{\mathrm{e}}\; {\mu }_{i}\; or\; {Adj\mu }_{i}\; {\mathrm{i}}{\mathrm{s}}\; {\mathrm{d}}{\mathrm{i}}{\mathrm{f}}{\mathrm{f}}{\mathrm{e}}{\mathrm{r}}{\mathrm{e}}{\mathrm{n}}{\mathrm{t}}\right),$ which is termed as:

The decision criteria for both performance measures are illustrated as follows: a test with the probability of the type Ⅰ error should be around α is declared the best test, while a test with a large power is deemed the best test. In order to give a quantitative definition of robustness (of significance level), we have to state the range of values of probability of type Ⅰ error for a given α value, for which the test would be considered robust. Bradley ^[72] suggested that a method could be regarded as robust to the violation of assumptions if the type Ⅰ error rate is within $\pm 0.5\alpha$. Bradley liberal criterion for robustness is $(0.5\alpha \le {\alpha }^{{\mathrm{*}}}\le 1.5\alpha)$. When $\alpha = 5{\mathrm{\%}}$, the estimated error rate outside the range $(2.5{\mathrm{\%}}, \; 7.5{\mathrm{\%}})$ is considered as conservative or liberal. Bradley's stringent criterion of robustness is $(0.9\alpha \le {\alpha }^{{\mathrm{*}}}\le 1.1\alpha)$.

Sullivan and D Agostino ^[73] reported a procedure as robust if the actual significance level does not exceed 10% of the nominal significance level $({\mathrm{e}}.{\mathrm{g}}., \; {\mathrm{f}}{\mathrm{o}}{\mathrm{r}}\; {\mathrm{\alpha }}\; = \; 0.05, \; {\mathrm{l}}{\mathrm{e}}{\mathrm{s}}{\mathrm{s}}\; {\mathrm{t}}{\mathrm{h}}{\mathrm{a}}{\mathrm{n}}\; {\mathrm{o}}{\mathrm{r}}\; {\mathrm{e}}{\mathrm{q}}{\mathrm{u}}{\mathrm{a}}{\mathrm{l}}\; {\mathrm{t}}{\mathrm{o}}\; \; \; \; 0.055)$. According to Guo and Luh ^[74], a method is robust if its observed significance level does not exceed 0.075 for the 5 percent nominal significance level. Zumbo and Coulmbo ^[75] expanded Bradley's robust criterion to identify three different levels of robustness. For $\alpha = 0.05$, the fairly stringent criterion is $(0.045, \; 0.055)$, the moderate criterion is $(0.04, \; 0.06)$, the liberal criterion is $(0.025, \; 0.075)$. Another criterion used by Vorapongsathorn et al. ^[76] is the Cochran limit, i.e., $(0.04, \; 0.06)$ for 5 percent nominal significance level. As there exists sampling error or some natural variation; therefore, to account for sampling error associated with estimated type Ⅰ error rates, we used Bradley's liberal criterion, to establish sampling error ranges around $\alpha$ in this study.

4.1.   Null case with homogeneity of variances

The percentage type Ⅰ error rates of both tests under the null case ${\delta }_{1}$ (no mean shift in any treatment) and variance homogeneity ${v}_{1}$(equal variance for all treatments) were reported in Table 6. The percentage type Ⅰ error that lies outside Bradley's liberal criterion range are tagged with the symbol ${\mathrm{"}}{\mathrm{*}}{\mathrm{"}}$. The findings of the current setup are listed below:

1) Under the balanced design, the ANOMC and ANOM tests have almost similar percentage type Ⅰ error rates except in heavy-tailed distributions such as exponential and Cauchy.

2) Similar findings are also observed under the unbalanced design, but the ANOMC test has an excessive percentage type Ⅰ error rate under all (normal or non-normal) environments when there is a direct pairing of correlations $\left({\rho }_{2}\right)$, while opposite results are observed in the case of the indirect pairing of the correlations $\left({\rho }_{3}\right)$.

3) Overall, the ANOMC test is not robust compared to the ANOM test when the response variable follows large heavy-tailed distributions.

4) Unequal correlations, either direct or indirect, may cause a change in the percentage type Ⅰ error rate from the specified $\alpha = 5{\mathrm{\%}}$.

5) Both tests reveal an approximately similar percentage of type Ⅰ error with the increase in the number of treatment levels and sample size.

4.2.   Null case with the heterogeneity of variances

As mentioned in Section 2.1, the ANOM test requires an assumption about the homogeneity of variances. Moreover, the ANOMC also works under the assumption that variances of conditional study variable are the same for each treatment group (cf. Section 2.2). In the ANOMC test, homogeneity may be categorized as; (i) the variances of $Y$ are equal for each treatment group (homogeneity), and (ii) the variances of $Y$ do not depend on the values of covariate $X$ (heteroscedasticity). In this study, we are concerned about the first condition of homogeneity, which significantly impacts the test performance under unbalanced design case. Therefore, we have introduced direct (i.e., ${v}_{2}^{+}, {v}_{3}^{+}$) and indirect (i.e., ${v}_{2}^{-}, {v}_{3}^{-}$) variance ratios to check the effect of heterogeneity on ANOMC and ANOM tests.

The impact of heterogeneity on ANOMC and ANOM tests with respect to different correlation pairs, sample sizes, distributional environments, and treatment levels $\left(t = 3, \; 4\; {\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; 5\right)$ have been investigated in this study and are reported in Tables 6–9. It is noted that type-Ⅰ error rates for ANOMC are slightly better than ANOM, but they exceed the nominal level at some stages (which is the effect of heterogeneity), although we have observed improvements at various levels of heterogeneity.

To be more specific, the prime findings of the effect of heterogeneity on the ANOMC and ANOM tests are listed below:

1) Under the balanced design, both tests (ANOMC and ANOM) are affected due to heterogeneity of variances, but the ANOMC test is less affected than the ANOM test in normal and non-normal environments except for exponential distribution.

2) Under the unequal sample sizes (unbalanced design), the same findings are still valid.

3) The ANOMC test has a higher type Ⅰ error rate (%) when large sample sizes are associated with more substantial variances, while the ANOM test has a higher type Ⅰ error rate (%) in the presence of an inverse relationship between sample sizes and variances.

The performance of both tests (ANOMC and ANOM) under heterogeneity is decreased with the increase in heterogeneity level. Meanwhile, when correlations are equal (${\rho }_{1}{\mathrm{a}}{\mathrm{n}}{\mathrm{d}}{\rho }_{4}$), the ANOMC test may produce relatively same percentage type Ⅰ error rates.

4.3.   Non-null cases under homogeneity of variances

For the null case, data has been sampled from a common population $\left({\mu }_{1} = {\mu }_{2} = \dots = {\mu }_{t}\right),$ and hence, any significance between treatment means attributed as sampling error and measured in terms of percentage type Ⅰ error $\left(\alpha \right)$. The non-null case consists of data that has been sampled from a population having at least one different mean, and the significance between the treatment means is measured in terms of the percentage power of the test $\left(1-\beta \right)$. In this study, nine different non-null cases (${\delta }_{2-10}$) are studied to examine the power of ANOMC and ANOM tests. Under the homogeneity of variances, the effect of several non-null cases on ANOMC and ANOM tests with respect to different correlation choices, distributional environments, and treatment levels $\left(t = 3, \; 4\; {\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; 5\right)$ are given in Table 10 for ${n}_{1}$, Table 11 for ${n}_{4}$ and Figures 1–4 for ${n}_{2}$, ${n}_{3}$, ${n}_{5}$ and ${n}_{6}$.

4.3.1.   Under the balanced design

The effect of several non-null cases on the ANOMC and ANOM tests with respect to distributional environments, treatment levels, and sample $\left({n}_{1}, {n}_{2}{\mathrm{a}}{\mathrm{n}}{\mathrm{d}}{n}_{3}\right)$ are given in Table 10 and Figures 1 and 2.

At the fixed sample size (n₁): The findings of the ANOMC and ANOM tests at the fixed sample size ${n}_{1}$ are reported in Table 10. At fixed correlations ${\rho }_{1}$, treatment level $\left(t = 3\right)$ and under Cauchy distribution, the findings of the non-null case ${\delta }_{2}$ reveals that the ANOMC test has 46.3% power as compared to 45.8% power of the ANOM test. However, the ANOMC and ANOM tests have 32.7% and 28.9% power for the non-null case ${\delta }_{4}$ under the exponential distribution. Further, under the double exponential distribution, the ANOMC and ANOM tests have 22.3% and 25.9% power for the non-null case ${\delta }_{6}$ and under the normal distribution; the ANOMC test has 18.0% power as compared to 21.9% power of the ANOM test for the non-null case ${\delta }_{8}$.

At the fixed sample size (n₂): The comparative analysis of the ANOMC and ANOM tests based on several non-null cases for the sample size choice ${n}_{2}$ are exhibited in Figure 1. At fixed $t = 4$, ${\rho }_{3},$ and ${\delta }_{3}$ under the normal distribution, the findings depict that the ANOMC test has 34.0% power as compared to 28.6% power of the ANOM test. The ANOMC and ANOM tests have 33.1% and 27.8% power for the non-null case ${\delta }_{5}$ under the double exponential distribution, while under the exponential distribution, the ANOMC and ANOM tests have 43.1% and 30.4% power for the non-null case ${\delta }_{7}$. Furthermore, for the non-null case ${\delta }_{9}$ under the Cauchy distribution, the ANOMC test has 47.2% power as compared to 36.1% power of the ANOM test.

At the fixed sample size (n₃): Several non-null cases for the ANOMC and ANOM tests at a fixed sample size ${n}_{3}$ are presented in Figure 2. When the treatment level $\left(t = 5\right)$ and correlations $\left({\rho }_{4}\right)$ are fixed than the findings of the non-null case ${\delta }_{2}$ shows that under double exponential distribution, the ANOMC test has 99.6% power as compared to 81.3% power of the ANOM test. Further, the ANOMC and ANOM tests have 98.4% and 67.4% power for the non-null case ${\delta }_{4}$ under the normal distribution. Under the exponential distribution, the ANOMC and ANOM tests have 99.7% and 82.9% power for the non-null case ${\delta }_{6}$ while for the non-null case ${\delta }_{10}$ under the Cauchy distribution, the ANOMC test has 98.1% power as compared to 78.7% power of the ANOM test.

4.3.2.   Unbalanced design

The comparative analysis of the ANOMC and ANOM tests with respect to several non-null cases, distributional environments, treatment levels, and sample sizes $\left({n}_{4}, {n}_{5}{\mathrm{a}}{\mathrm{n}}{\mathrm{d}}{n}_{6}\right)$ are given in Table 11 and Figures 3 and 4.

At the fixed sample size (n₄): The effect of several non-null cases on the performance of ANOMC and ANOM tests with respect to distributional environments and treatment levels are reported in Table 11. On the fixed treatment level $\left(t = 5\right)$, correlations $\left({\rho }_{4}\right)$, non-null case $\left({\delta }_{2}\right),$ and under double exponential distribution, the ANOMC test has 95.4% power as compared to 71.2% power of the ANOM test. Under the normal distribution, the ANOMC and ANOM tests have 69.7% and 42.8% power for the non-null case ${\delta }_{4}$ and under the exponential distribution, the ANOMC and ANOM tests have 85.5% and 57.2% power for the non-null case ${\delta }_{6}$. Further, for the non-null case ${\delta }_{10}$ and under the Cauchy distribution, the ANOMC test has 65.1% power as compared to 38.2% power of the ANOM test.

At the fixed sample size (n₅): The performance analysis of the ANOMC and ANOM test at a fixed sample size ${n}_{5}$ are presented in Figure 3. For the fixed treatment level $\left(t = 4\right)$ having the fixed correlations $\left({\rho }_{3}\right)$, the findings of the non-null case ${\delta }_{3}$ depicts that under the normal distribution, the ANOMC test has 49.7% power as compared to 41.3% power of the ANOM test. The ANOMC and ANOM tests have 41.5% and 33.9% power for the non-null case ${\delta }_{5}$ under the double exponential distribution. Further, under the exponential distribution, the ANOMC and ANOM tests have 81.5% and 68.5% power for the non-null case ${\delta }_{7}$ while for the non-null case ${\delta }_{9}$ under the Cauchy distribution, the ANOMC test has 82.0% power as compared to 71.9% power of the ANOM test.

At the fixed sample size (n₆): Several non-null cases for the ANOMC and ANOM tests at a fixed sample size ${n}_{6}$ are presented in Figure 4. At fixed equal correlations $\left({\rho }_{2}\right)$ having the treatment level $\left(t = 3\right)$, the findings of the non-null case ${\delta }_{2}$ reveal that under the Cauchy distribution, the ANOMC test has 94.7% power as compared to 85.8% power of the ANOM test. Further, under the exponential distribution, the ANOMC and ANOM tests have 72.7% and 53.6% power for the non-null case ${\delta }_{4}$ and under the double exponential distribution, the ANOMC and ANOM tests have 59.9% and 43.6% power for the non-null case ${\delta }_{6}$ while for the non-null case ${\delta }_{8}$ under the normal distribution, ANOMC test has 47.5% power as compared to 35.0% power of the ANOM test.

Overall, the performance of both tests (ANOMC and ANOM) increases with the increase in sample sizes, but the ANOMC test performs relatively better than the ANOM test. Moreover, the performance of both tests (ANOMC and ANOM) also increases due to the increase in correlation between the study variable and the concomitant variable. Furthermore, when sample sizes are unequal, the performance of both tests is affected by the inverse relationship between sample sizes and the size of shifted means.

5.   Experimental examples

In this section, two illustrative examples from different experimental situations are discussed to compare the performance of the proposed method ANOMC and the ANOM test.

5.1.   An illustrative example of the balanced case

For equal sample sizes (balanced case), the ANOM and ANOMC methods are implemented on the mechanical manufacturing problem dataset. Electrical discharge machining (EDM) is a frequently used method in the manufacturing industry. Dutta et al. ^[54] described an experimental study to investigate the effects of EDM parameters (i.e., pulse current, pulse-on-time, and pulse-off-time) on machining time and surface roughness for machining Inconel 800. The experimental work was carried out on the Electronic4-axis CNC sprint cut wire electrical discharge machine. A negatively polarized brass wire of diameter 0.25mm with a tensile strength of 500N/mm was used as an electrode. Deionized water was used as the dielectric fluid. Samples of size 25mm×25mm×5mm were cut on the machine, and the machining time (min) and surface roughness $\left(\mu m\right)$ with the pulse current $\left(amp\right)$ are reported therein. In this example, we are considering surface roughness as the study variable, machining time as a concomitant variable and both variables are reported with several levels of pulse current $\left(amp\right)$. It is noted that surface roughness has a linear relation with machining time based on ten experiments, excluding the first and seventh experiments. Therefore, 10 observations $\left(n = 10\right)$ with respect to the three levels $\left(t = 3\right)$ of pulse current (i.e., $210\mu m, \; 220\mu m$ and $230\mu m)$ are used to implement ANOM and ANOMC methods.

The ANOM method is applied to the observations of surface roughness without incorporating the concomitant variable machining time, and the results are plotted in Figure 5. The overall average $\left(\overline{{Y}_{..}}\right)$ and mean square error $\left(MSE\right)$ of surface roughness are calculated as 3.072 and 0.318, respectively. Using Table 2, the critical value (i.e., $h\left(\alpha, t, N-t\right) = 2.51$) is fixed against the level of significance $\alpha = 5{\mathrm{\%}}$. Further, the individual means of pulse current levels are plotted against the decision interval (i.e., $UDL = 2.706\; {\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; LDL = 3.437$) gives evidence of not rejecting the null hypothesis, i.e., no individual average differs from the overall average. Hence, the power of all pulse current levels to detect the difference in surface roughness is likely small.

Further, the ANOMC method is applied for testing the means of surface roughness without ignoring the effect of machining time. The surface roughness and machining time correlations are calculated as -0.964433, -0.9618967, and -0.8263165 with respect to pulse current levels. Moreover, the regression means are calculated for each pulse current level, and the overall regression average $\left({\overline{M}}_{..}\right),$ and mean square error of regression mean estimator $\left({MSE}_{{\overline{M}}_{..}}\right)$ are estimated as 3.073646 and 0.07143622, respectively. Using Table 2, the critical value (i.e., ${h}^{{\mathrm{*}}}(\alpha, t, n) = 8.6$) is fixed against the level of significance $\alpha = 5{\mathrm{\%}}$. Further, the graphical representation of the ANOMC under balanced design is also plotted in Figure 5, where individual regression means are plotted against the decision interval (i.e., $UDL = 2.915021\; {\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; LDL = 3.232271$). The findings reveal that all individual regression averages differ from the overall regression average, which is evidence of the significance of surface roughness with respect to pulse current levels.

5.2.   An illustrative example of the unbalanced case

For the unbalanced design (unequal sample sizes), dataset related to medical science is used, where the effect on participants and partner libidos are reported with respect to three Viagra dosages. The complete data set having three treatments (Viagra dosages) with different sample sizes of participants and partner libido is reported in Table 11.1 on page no. 400 ^[77]. The first Viagra dosage (Placebo) has 9 samples of participant and partner libidos, while other Viagra dosages (Low dose and High dose) have 8 and 13 samples of participant and partner libidos, respectively.

In this example, we used the number of participant libido as a study variable $\left(Y\right),$ while the number of partner libido is used as a covariate $\left(X\right)$. The graphical layout of the ANOM method is presented in Figure 6. The ANOM test is applied to the participant's libido observations by ignoring the partners' libido effect. Using Table 2, the critical value $m\left(\alpha, t, N-t\right) = 2.53$ is fixed against the $5{\mathrm{\%}}$ level of significance. The overall average $\left(\overline{{Y}_{..}}\right),$ and mean square error $\left(MSE\right)$ of participant libidos are calculated as 4.366667 and 3.486032, respectively. The individual means of Viagra dosage are plotted against the decision intervals, i.e., $UDL = \left(5.684055;5.796851;5.352897\right){\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; LDL = (3.049278;\; 2.936483;3.380436)$), which reveals that all individual averages are statistically insignificant. Hence, all Viagra dosages have a similar effect on participant's libido.

The ANOMC method is applied for testing the means of participants' libido without ignoring the effect of partners' libido. The correlations of participant and partner libidos are calculated as 0.8829347, 0.9718268, and -0.1688756 with respect to Viagra dosages. Using Table 2, the critical value ${h}^{{\mathrm{*}}}(\alpha, t, n) = 8.63$ is fixed against the $5{\mathrm{\%}}$ level of significance. Moreover, the regression means are calculated for each Viagra dosage, and the overall regression average $\left({\overline{M}}_{..}\right)$, and mean square error of regression mean estimator $\left({MSE}_{{\overline{M}}_{..}}\right)$ are estimated as 4.641178 and 0.3935578, respectively. The graphical representation of the ANOMC method under balanced design is also plotted in Figure 6, where individual regression means are plotted against the decision interval, i.e., $UDL = \left(5.588390;5.669491;5.350285\right){\mathrm{a}}{\mathrm{n}}{\mathrm{d}}\; LDL(3.693966;3.612866;3.932072)$. The findings reveal that only individual regression averages related to the placebo drug are different from the overall regression average, which is evidence that a placebo drug has a different effect on the participant's libido.

6.   Conclusions

ANOVA is the most commonly used technique to compare the treatment means. An alternative technique to ANOVA is ANOM, a graphical test used to test whether the treatment means differ from the grand mean. ANOVA requires multiple comparison tests to identify the significantly different treatments; however, ANOM does not require any additional test for such identification. This study proposed a new covariate based ANOM method, namely ANOMC, for the analysis of means. It is used for testing the significance of means from the grand mean by accommodating the effect of a covariate. The proposed procedure works under several assumptions, such as normality, linearity, and homogeneity. The effect of these assumptions, sample sizes (equal or unequal), treatments, and hypotheses (null and non-null) on ANOM and ANOMC tests are compared in terms of percentage type Ⅰ error and percentage power of the test.

The findings of the study revealed that in the case of homogeneity of variances with the null case, the ANOMC test is not as robust as the ANOM test when the response variable follows a conditionally large heavy-tailed distribution (e.g., exponential distribution). It is observed that under unequal correlations, pairing (direct or indirect) of correlations may lead to a change in the percentage type Ⅰ error rate from pre-specified $\alpha = 5{\mathrm{\%}}$. Moreover, both tests have approximately similar findings with the increase in treatment level and sample size. In the presence of heterogeneity of variances with the null case, both tests are affected, but the ANOMC test is less affected as compared to the ANOM test in a balanced design under normal and non-normal environments except for exponential distribution. In an unbalanced design, the ANOMC test is affected when large sample sizes are associated with more substantial variances, while the ANOM test has a higher type Ⅰ error rate in the presence of an inverse relationship between sample sizes and variances. It is also noted that the ANOMC test has relatively the same percentage type Ⅰ error rate for equal correlations.

As expected, the power values of the ANOM test and ANOMC test change with respect to effect size $\left({\mathrm{\delta }}\right)$, sample size, treatment level, and distribution environment. The performance of both tests improves with the increase in sample sizes, but the ANOMC test performs relatively better than the ANOM test. The performance of both tests (ANOMC and ANOM) also increases due to the increase in correlation between the study variable and the concomitant variable. Moreover, when sample sizes are unequal, the performance of both tests is affected by the inverse relationship between sample sizes and the size of shifted means. This study is designed under a limited number of treatments, choices of sample size (equal and unequal), correlations, and distributional environments, which may be extended in the future. Moreover, the proposal may also be expanded using the robust regression estimators to achieve a robust version of the ANOMC method. In the current study, we have used Monte Carlo simulations. However, developing other tests (integral approach or Markov chain method) for ANOMC to construct LDL and UDL is a potential direction for future research.

Acknowledgment

This research work is based upon work supported by King Fahd University of Petroleum & Minerals. Authors at KFUPM acknowledge the Interdisciplinary Research Center for Smart Mobility and Logistics for the support received under Grant no. INML2205.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

R code for example about EDM problem

y = c (2.05, 2.43, 2.79, 2.85, 3.14, 2.12, 2.85, 3.35, 3.19, 3.38, 2.1, 2.79, 3.24, 3.43, 3.59,

2.35, 3.34, 3.4, 3.44, 4.02, 2.56, 2.81, 3.17, 3.44, 4.16, 2.31, 2.74, 3.49, 3.65, 3.97)

x = c (31.67, 19.25, 14.5, 12.15, 10, 26.6, 14, 9.47, 8.88, 7.87, 30.27, 18.3, 13.6, 11.87, 10.1,

25.4, 13.4, 8.93, 8.07, 7.67, 29.67, 15.28, 13.4, 11.58, 9.07, 25, 12.3, 8.4, 7.83, 6.13)

t = 3; n = 10; c1 = 2.51; c2 = 8.6

ymat = matrix (y, n, t); xmat = matrix (x, n, t)

par (mfrow = c (1, 2))

par (mar = c (2, 4.2, 2, 3), cex = 1)

################################# ANOM Test#################################

ybari = apply (ymat, 2, mean); s2i = apply (ymat, 2, var)

ybar = mean (ybari); mse = mean (s2i)

ldl = ybar-(c1*sqrt (mse)*sqrt ((t-1)/(n*t)))

udl = ybar + (c1*sqrt(mse)*sqrt ((t-1)/(n*t)))

plot (ybari, ylim = c (2.5, 3.5), col = "blue", pch = 20, cex = 1.5, xlab = "Treatments",

ylab = expression (bar (Y~scriptstyle(i.))), xaxt = "n", main = "ANOM")

axis (labels = list ("210 (amp)", "220 (amp)", "230 (amp)"), side = 1, at = c (1, 2, 3), cex.axis = 0.75)

g = c ()

for (i in 1:t) {

if(ybari[i] > udl|ybari[i] < ldl){g[i] = ybari[i]} else {g[i] = NA}

}

points (g, cex = 1.5, col = "red", pch = 20, lty = 2)

abline (h = ldl, v = NULL, col = "green4", lty = 5)

abline (h = ybar, v = NULL, col = "green4", lty = 1)

abline (h = udl, v = NULL, col = "green4", lty = 5)

segments (1, ybar, 1, ybari[1], col = "blue")

segments (2, ybar, 2, ybari[2], col = "blue")

segments (3, ybar, 3, ybari[3], col = "blue")

mtext (("LDL"), side = 4, line = 1, at = ldl, cex = 0.75, col = "Green4")

mtext (expression(paste(bar(Y[..]))), side = 4, line = 1, at = ybar, cex = 0.75, col = "Green4")

mtext(("UDL"), side = 4, line = 1, at = udl, cex = 0.75, col = "Green4")

############################ANOMC Test #####################################

sdy = apply (ymat, 2, sd); sdx = apply (xmat, 2, sd)

ybari = apply (ymat, 2, mean); ybar = mean (ybari)

xbari = apply (xmat, 2, mean); xbar = mean (xbari)

rr = beta = mm = k = smm = double ()

for (i in 1:t)

{

rr[i] = cor(ymat[, i], xmat[, i])

beta[i] = rr[i]*(sdy[i]/sdx[i])

mm[i] = ybari[i]-beta[i]*(xbar-xbari[i])

k[i] = sqrt((1-(rr[i]^2))*(1+(1/(n-3))))

smm[i] = (k[i]*sdy[i])/sqrt(n)

}

mbar = mean (mm); mse = mean (smm)

ldl = mbar-(c2*mse*sqrt((t-1)/(n*t)))

udl = mbar+(c2*mse*sqrt((t-1)/(n*t)))

plot (mm, ylim = c (2.5, 3.5), col="blue", pch = 20, cex = 1.5, xlab="Treatments",

ylab = expression (bar(M~scriptstyle(i.))), xaxt = "n", main="ANOMC")

axis (labels = list("210(amp)", "220(amp)", "230(amp)"), side = 1, at = c(1, 2, 3), cex.axis = 0.75)

g = c ()

for (i in 1:t){

if(mm[i] > udl|mm[i] < ldl){g[i] = mm[i]} else {g[i] = NA}

}

points (g, cex = 1.5, col="red", pch = 20, lty = 2)

abline (h = ldl, v = NULL, col = "green4", lty = 5)

abline (h = mbar, v = NULL, col = "green4", lty = 1)

abline (h = udl, v = NULL, col = "green4", lty = 5)

segments (1, mbar, 1, mm[1], col="blue")

segments (2, mbar, 2, mm[2], col="blue")

segments (3, mbar, 3, mm[3], col="blue")

mtext (("LDL"), side = 4, line = 1, at = ldl, cex = 0.75, col = "Green4")

mtext (expression(paste(bar(M[..]))), side = 4, line = 1, at = mbar, cex = 0.75, col = "Green4")

mtext (("UDL"), side = 4, line = 1, at = udl, cex = 0.75, col = "Green4")

R code for example about viagra dosages problem

p0 = c(3, 2, 5, 2, 2, 2, 7, 2, 4); p1 = c(4, 1, 5, 1, 2, 2, 7, 4, 5)

l0 = c(7, 5, 3, 4, 4, 7, 5, 4); l1 = c(5, 3, 1, 2, 2, 6, 4, 2)

h0 = c(9, 2, 6, 3, 4, 4, 4, 6, 4, 6, 2, 8, 5); h1 = c(1, 3, 5, 4, 3, 3, 2, 0, 1, 3, 0, 1, 0)

n1 = length(p0); n2 = length(l0); n3 = length(h0);

t = 3; n = c(n1, n2, n3); N = sum(n); c1 = 2.53; c2 = 8.63

par (mfrow = c(1, 2))

par (mar = c(2, 4.2, 2, 3), cex = 1)

######################### ANOM Test #####################################

ybari = c(mean(p0), mean(l0), mean(h0))

s2i = c(var(p0), var(l0), var(h0))

ybar = ((n[1]*ybari[1])+(n[2]*ybari[2])+(n[3]*ybari[3]))/N

mse = (((n[1]-1)*s2i[1])+((n[2]-1)*s2i[2])+((n[3]-1)*s2i[3]))/(N-t)

ldl = udl = double()

for (k in 1:t)

{

ldl[k] = ybar-(c1*sqrt(mse)*sqrt((N-n[k])/(N*n[k])))

udl[k] = ybar+(c1*sqrt(mse)*sqrt((N-n[k])/(N*n[k])))

}

plot (ybari, ylim = c(2.5, 6.5), col = "blue", pch = 20, cex = 1.5, xlab = "Treatments",

ylab = expression(bar(Y~scriptstyle(i.))), xaxt = "n", main="ANOM")

axis (labels = list("Placebo", "Low dose", "High dose"), side = 1, at = c(1, 2, 3), cex.axis = 0.75)

g = c()

for (i in 1:t){

if(ybari[i] > udl[i]|ybari[i] < ldl[i]){g[i] = ybari[i]} else {g[i] = NA}

}

points (g, cex = 1.5, col = "red", pch = 20, lty = 2)

abline (h = ybar, v = NULL, col = "green4", lty = 1)

segments (0, udl[1], 1.5, udl[1], col = "green4", lty = 5)

segments (1.5, udl[2], 2.5, udl[2], col = "green4", lty = 5)

segments (2.5, udl[3], 3, udl[3], col = "green4", lty = 5)

segments (0, ldl[1], 1.5, ldl[1], col = "green4", lty = 5)

segments (1.5, ldl[2], 2.5, ldl[2], col = "green4", lty = 5)

segments (2.5, ldl[3], 3, ldl[3], col = "green4", lty = 5)

segments (1, ybar, 1, ybari[1], col = "blue")

segments (2, ybar, 2, ybari[2], col = "blue")

segments (3, ybar, 3, ybari[3], col = "blue")

mtext (("LDL"), side = 4, line = 1, at = mean(ldl), cex = 0.75, col = "Green4")

mtext (expression(paste(bar(Y[..]))), side = 4, line = 1, at = ybar, cex = 0.75, col = "Green4")

mtext (("UDL"), side = 4, line = 1, at = mean(udl), cex = 0.75, col = "Green4")

################################ANOMC Test ##################################

ybari = c(mean(p0), mean(l0), mean(h0));

xbari = c(mean(p1), mean(l1), mean(h1)); xbar = mean(xbari)

sdy = c(sd(p0), sd(l0), sd(h0))

sdx = c(sd(p1), sd(l1), sd(h1))

rr1 = cor(p0, p1); rr2 = cor(l0, l1); rr3 = cor(h0, h1); rr = c(rr1, rr2, rr3)

beta = mm = k = smm = double ()

for (i in 1:t)

{

beta[i] = rr[i]*(sdy[i]/sdx[i])

mm[i] = ybari[i]-beta[i]*(xbar-xbari[i])

k[i] = sqrt((1-(rr[i]^2))*(1+(1/(n[i]-3))))

smm[i] = (k[i]*sdy[i])/sqrt(n[i])

}

mbar = ((n[1]*mm[1])+(n[2]*mm[2])+(n[3]*mm[3]))/N

mse = (((n[1]-1)*smm[1])+((n[2]-1)*smm[2])+((n[3]-1)*smm[3]))/(N-t)

ldl = udl = double()

for (k in 1:t)

{

ldl[k] = mbar-(c2*mse*sqrt((N-n[k])/(N*n[k])))

udl[k] = mbar+(c2*mse*sqrt((N-n[k])/(N*n[k])))

}

plot (mm, ylim = c(3, 6), col = "blue", pch = 20, cex = 1.5, xlab = "Treatments",

ylab = expression(bar(M~scriptstyle(i.))), xaxt = "n", main = "ANOMC")

axis (labels = list("Placebo", "Low dose", "High dose"), side = 1, at = c(1, 2, 3), cex.axis = 0.75)

g = c ()

for (i in 1:t){

if (mm[i] > udl[i]|mm[i] < ldl[i]){g[i] = mm[i]} else {g[i] = NA}}

points (g, cex = 1.5, col = "red", pch = 20, lty = 2)

abline (h = mbar, v = NULL, col = "green4", lty = 1)

segments (0, udl[1], 1.5, udl[1], col = "green4", lty = 5)

segments (1.5, udl[2], 2.5, udl[2], col = "green4", lty = 5)

segments (2.5, udl[3], 3, udl[3], col = "green4", lty = 5)

segments (0, ldl[1], 1.5, ldl[1], col = "green4", lty = 5)

segments (1.5, ldl[2], 2.5, ldl[2], col = "green4", lty = 5)

segments (2.5, ldl[3], 3, ldl[3], col = "green4", lty = 5)

segments (1, mbar, 1, mm[1], col = "blue")

segments (2, mbar, 2, mm[2], col = "blue")

segments (3, mbar, 3, mm[3], col = "blue")

mtext (("LDL"), side = 4, line = 1, at = mean(ldl), cex = 0.75, col = "Green4")

mtext (expression(paste(bar(M[..]))), side = 4, line = 1, at = mbar, cex = 0.75, col = "Green4")

mtext (("UDL"), side = 4, line = 1, at = mean(udl), cex = 0.75, col = "Green4")