# Detecting Changes under Multivariate Normal Distributions via the Generalized Inference

Authors: Weiyan Mu; Xin Wang; Xi Wu; Shifeng Xiong
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

It is commonly encountered in many fields to detect whether a change occurs on a population after a special process. Based on observations for describing the population before and after the process, we formulate this problem as two statistical hypotheses testing problems within a framework of multivariate statistical analysis and then propose a generalized inference approach to solve them. The corresponding generalizedp values and their calculation details are provided. The proposed method is also extended to multiple testing problems. Simulation studies show that the proposed p values have satisfactory frequentist performance. We illustrate our methods with a real application in manufacturing of bearings that are used in medical devices.


---

## Body

## 1. Introduction

In many practical problems, we need to analyze uncertainty of a population induced by an operation. Specifically, letX∈ℝp be the p indices that describe the population. For n independent individuals from the population, let X1,…,Xn denote their values of X. After some operation acts on these individuals, their values become Y1,…,Yn, respectively. These Y1,…,Yn can be viewed as an independent sample from Y∈ℝp, the p indices of the population after the operation. We are interested in whether the operation has a significant influence on the population, i.e., there exists a significant change from X to Y.There are a number of examples of the above problem in real applications. In medicine, the population can be the whole patients, and their indicesX become Y after they take a kind of drugs. In manufacturing, the population can be the whole products, and their indices X become Y after a certain process. We want to know whether the change from X to Y is statistically significant.Note thatp is usually larger than 1, i.e., we need to consider more than one indices simultaneously. This paper studies the problem of testing changes within a framework of multivariate statistical analysis. Let Z=Y−X, and assume Z∼ a multivariate normal distribution Npμ,Σ with unknown mean μ and covariance matrix Σ, where Σ is positive definite. The value of θ=EY−X2=EZ2=μ2+trΣ can be used to quantify the uncertainty of the change from X to Y, where tr denotes the trace of a matrix, and we let it be the parameter of interest. Let the independent sample of Z be Z1,…,Zn, where Zi=Yi−Xi for i=1,…,n. With a prespecified θ0>0, the problem of testing changes is formulated as the following statistical hypothesis testing problems:(1)H0:θ⩽θ0,↔H1:θ>θ0,(2)H0:θ⩾θ0,↔H1:θ<θ0.where θ0 is prespecified by experienced practitioners and can be viewed as a threshold value for justifying whether a change significantly occurs from X to Y. In other words, we say that a change significantly occurs if and only if θ is greater than θ0. When we need strong evidence to support that change occurs, we should consider testing (1). This can be used in testing whether a drug is effective. When we need strong evidence to support that there is no change, we should consider testing (2). This can be used in manufacturing where we expect that some process does not influence the product.It should be pointed out that the above statistical problems are completely different from the conventional two-sample comparison test sinceX is not independent of Y. They are indices of the same population in two phases. Problems (1) and (2) are also different from mean tests or covariance tests of multivariate normal distributions, which can be derived based on the Wilks statistics (likelihood ratio statistics) [1], since the parameter θ involves both the mean and a part of the covariance matrix. To the best of our knowledge, there is no work on the above statistical hypothesis testing problems for detecting changes. Existing research studies on change detection in the literature [2–4] relate to different statistical models or focuses. Note that the distributions of the likelihood ratio statistics under the null one-sided hypotheses like (1) and (2) are very complicated [5]. In this paper, we use the generalized inference approach [6, 7] to provide generalized p-values of testing (1) and (2). The concept of the generalized p value was introduced by Tsui and Weerahandi [6]; this idea was extended to interval estimation by Weerahandi [7]. Since then the generalized inference approach has been successfully applied to many complex inferential problems, many of which appear in industry and medical science, see Mu et al. [8], Chen and Ye [9], and Feng and Tian [10] among many others. For multivariate normal distributions, Gamage et al. [11] and Park [12] proposed the generalized inference methods for specified problems, which are different from the above testing problems in (1) and (2). The generalized inference approach is proven to be asymptotically correct under mild conditions [13–15]. Simulation results in the literature indicate that this approach usually possesses good frequentist properties. Another appealing feature of the generalized inference approach is its simplicity under some parametric models, especially normal models. We do not need to derive complex (asymptotic) distributions of test statistics such as the likelihood ratio statistics for one-sided hypotheses.This paper is organized as follows. Section2 presents the generalized inference approach for testing (1) and (2). Section 3 extends the proposed method to multiple testing problems. In Section 4, we make a simulation study for Type-I error and power performance of the proposed methods. We apply our method to a dataset in manufacturing in Section 5. We conclude the paper with some concluding remarks in Section 6.


## 2. Methodology

Based on the sampleZ1,…,Zn independently and identically distributed from Npμ,Σ, we first construct the generalized pivotal quantities of the unknown parameters μ,Σ, which can induce the generalized pivotal quantity of θ. Then, the generalized p values for testing (1) and (2) can be given by the generalized pivotal quantity of θ.Let(3)Z¯=1n∑i=1nZi,V=∑i=1nZi−Z¯Zi−Z¯′.It is known thatV∼Wpn−1,Σ, the Wishart distribution with n−1 degrees of freedom and scale parameter matrix Σ. Let T be the Cholesky factor of V. That is, V=TT′, where T is a lower triangular matrix with positive diagonal elements. Such a decomposition is unique for any positive definite matrix [16]. Let Θ be the Cholesky factor of Σ with positive diagonal elements. Write(4)U=Θ−1T,and we see that U is a lower triangular matrix with positive diagonal elements and UU′=Θ−1TT′Θ′−1∼Wpn−1,Ip, which indicates that U is the unique Cholesky factor of Wpn−1,Ip. On the contrary, we write(5)nΘ−1Z¯−μ=E,and E∼Np0,Ip is independent of U. Based on the construction method of generalized pivotal quantities introduced by [13], we can obtain from (4) and (5) that the generalized pivotal quantities of μ and Θ are(6)Rμ=Z¯−TU−1En,RΘ=TU−1.Therefore, the generalized pivotal quantity ofθ is(7)Rθ=Rμ2+trRΘRΘ′=Z¯−TU−1En2+trTU−1U′−1T′.By (7), the generalized p-values for testing (1) and (2) are given by(8)p=PRθ⩽θ0,(9)p=PRθ⩾θ0,respectively.Note thatRθ can be viewed as a distribution estimator [17] or a generalized bootstrap variable [15] for θ. Some experience shows that Rθ overestimates θ in finite-sample cases, and this will lead to biases in testing one-side hypotheses such as (1) and (2). To overcome the difficulty, here we use Xiong [15] generalized bootstrap variable method to present new p values for (1) and (2). Let(10)θ^=Z¯2+trVnbe the maximum likelihood estimator of θ. By Xiong [15] method, Rθ−θ^ can be used to approximate the distribution of θ^−θ0, and this leads to p values for testing (1) and (2):(11)p=PRθ⩾2θ^−θ0,(12)p=PRθ⩽2θ^−θ0,respectively. We can use some combinations of the above results to give p values. Algorithm 1 provides implemental details for computing the p value in (11), and those for computing (8), (9), and (12) can be obtained similarly.Algorithm 1: Computing the p value (11).
Inputs: the sampleZ1,…,Zn, the threshold value θ0, and the Monte Carlo sample size MSteps:(1)
ComputeZ¯ and V in (3) and θ^ in (10).(2)
Compute the Cholesky factorT of V with positive diagonal elements.(3)
Form=1,…,M:GenerateE∼Np0,Ip independent of W∼Wpn−1,In.Compute the Cholesky factorU of W with positive diagonal elements.By (7), compute the mth observation Rθ,m.Endm loop.(4)
Computep=M−1∑m=1MVRθ,m⩾2θ^−θ0, where I is the indicator function.


## 3. Extension to Multistage Detections

We sometimes need to detect the changes after a multistage process. LetX1 denote the initial indices of the population. After the kth stage, Xk becomes Xk+1 for k=1,…,K. Denote θk=EXk+1−Xk2=μk2+trΣk, where Xk+1−Xk∼Npμk,Σk with unknown mean μk and unknown positive definite covariance matrix Σk for k=1,…,K. We want to know which ones among the K changes, Xk+1−Xkk=1,…,K, are statistically significant. Then, problems (1) and (2) should be extended to the following two multiple testing problems:(13)H0k:θk⩽θ0k,↔H1k:θk>θ0k,k=1,…,K,and(14)H0k:θk⩾θ0k,↔H1k:θk<θ0k,k=1,…,K,respectively. Here, the threshold values θ0k,k=1,…,K, are prespecified.LetZk=Xk+1−Xk. Based on the independent sample Zi1,…,ZiKi=1,…,n, we provide a Bonferroni inequality method for (13) and (14). The Bonferroni inequality is commonly used in multiple comparison problems [18]. Algorithm 2 presents implemental details of our method for (13), and those for (14) can be obtained similarly. It can be seen that the Bonferroni inequality method is simple to implement. It does not depend on the correlations between Zk and can control the family-wise error rate ⩽α [18].Algorithm 2: The Bonferroni inequality method for (13).
Inputs:The sampleZi1,…,ZiKi=1,…,n, the threshold values θ0kk=1,…,K, and the nominal level α∈0,1.Steps:(1)
Fork=1,…,K:Use sampleZiki=1,…,n to compute the p-value pk for testingH0k:θk⩾θ0k,↔H1k:θk<θ0kby the generalized inference approach in Section2.Endk loop.(2)
Reject the null hypothesesH0k if pk<α/K.Inputs:The sampleZi1,…,ZiKi=1,…,n, the threshold values θ0kk=1,…,K, and the nominal level α∈0,1.Steps:(1)
Fork=1,…,K:Use sampleZiki=1,…,n to compute the p-value pk for testingH0k:θk⩾θ0k,↔H1k:θk<θ0kby the generalized inference approach in Section2.Endk loop.(2)
Reject the null hypothesesH0k if pk<α/K.


## 4. Simulation Study

We now illustrate the proposed procedure with some simulated data. We consider the following cases:(15)CaseI:μ=11,Σ=10.50.51.5,CaseII:μ=0.5−1,Σ=1−0.5−0.51,CaseIII:μ=0000,Σ=10.20.20.10.21−0.50.20.2−0.52−0.20.10.2−0.21,CaseIV:μ=−0.5010,Σ=10.20.20.10.2100.20.201−0.20.10.2−0.21.For each case, we use Algorithm1 to implement our tests when θ0=μ2+trΣ. The combinations of n,p are given in Table 1. The Monte Carlo sample size M in Algorithm 1 is set as 10,000. The empirical Type-I errors over 10,000 repetitions are shown in Table 1. We can see that most errors are close to or less than the nominal level α. As n increases, the errors become closer to α. These findings indicate that the proposed testing methods are effective.Table 1
Type-I errorsα=0.05.
Case (I)Case (II)Case (III)Case (IV)Testing (1), p=2Testing (2), p=2Testing (1), p=4Testing (2), p=4n=20n=30n=40n=50n=20n=30n=40n=50n=40n=60n=80n=100n=40n=60n=80n=1000.0610.0540.0560.0610.0270.0470.0470.0530.0980.0680.0580.0490.0060.0180.0250.0340.0730.0530.0540.0460.0340.0320.0450.0460.0660.0500.0380.0390.0160.0240.0380.045In addition, we compute the powers of our test for (1) under Case (I). First, we let n increase and compute the powers at μ=1.51 and Σ=10.60.62. The results are shown in the left side of Figure 1. Then, we fix n=30, and compute the powers at μ=1+s1+s and Σ=1+s0.50.51.5+s with increasing s. The results are displayed in the right side of Figure 1. It can be seen that the power performance of our method is reasonable: the power tends to one fast as n or s increases.Figure 1
Powers under case (I)α=0.05.


## 5. Case Study

A bearing is the basic element of any rotating mechanism [19]. For a certain type of bearings that are used in medical devices, we investigate the manufacturing process of their retainers. Figure 2 shows a type of bearing retainers. There is a postprocessing in the entire manufacturing process, including two cleaning stages and one heat treatment stage. The objective is to assess whether the inner and outer diameters of the retainers significantly change after each stage in the postprocessing.Figure 2
A type of bearing retainers.In this study, we obtain 34 sample retainers and measure their inner and outer diameters before and after the postprocessing. The original data cannot be open due to a confidentiality agreement, and thus, they were slightly modified via a linear transformation. We use the modified data to conduct our analysis. For theith retainer, let Xik=Xi1k,Xi2k′ and Xik+1=Xi1k+1,Xi2k+1′ denote its diameters before and after stage k=1,2,3, respectively, where the first and second components of Xik represent the inner and outer diameters, respectively. The scatter plots of Xik and Zik=Xik+1−Xik are presented in Figure 3 for k=1,2,3. Based on the data, we test the multiple hypotheses (14) using the Bonferroni inequality method in Algorithm 2. The values of θ0 in the three stages are, respectively, 6×10−5, 6×10−5, and 6×10−4, which are suggested by engineers. With the Monte Carlo sample size M=106, we obtain the generalized p values p=0.0431, p=0.0022, and p=7.02×10−4, respectively. Given the nominal level α=0.05, we should reject H02 and H03, and this indicates that the inner and outer diameters do not significantly change after the second cleaning stage and the heat treatment stage.Figure 3
Scatter plots (unit: mm). For eachk=1,2,3, the white circles, asterisks, and black dots represent Xik, Xik+1, and Zik, respectively.


## 6. Concluding Remarks

In this paper, we have proposed a multivariate statistical approach to detect whether a change occurs after some process. The parameter of interest in our approach describes both changes in mean and variance terms and can be acted as an index to detect the change of the population. The generalizedp values for testing whether the parameter of interest exceeds a threshold value have been presented, and we have shown them to have good performance in our simulations. We also apply our methods to analyze products of bearing retainers and conclude that they do not significantly change after two stages in the postprocessing. A future topic is to extend our methods to other models when the normal assumption does not match the data. For such cases, it may be difficult to use the generalized inference approach. We will consider to combine this approach and other popular inferential methods such as the bootstrap method [20, 21]. Another possible future direction is to study the false discovery rate control method [22] when the number of stages in Section 3 is large.



---

*Source: 5526717-2021-05-25.xml*