# Tolerance Intervals in a Heteroscedastic Linear Regression Context with Applications to Aerospace Equipment Surveillance

Authors: Janet Myhre; Daniel R. Jeske; Michael Rennie; Yingtao Bi
Journal: International Journal of Quality, Statistics, and Reliability
 (2009)

---

## Abstract

A heteroscedastic linear regression model is developed from plausible assumptions that describe the time evolution of performance metrics for equipment. The inherited motivation for the related weighted least squares analysis of the model is an essential and attractive selling point to engineers with interest in equipment surveillance methodologies. A simple test for the significance of the heteroscedasticity suggested by a data set is derived and a simulation study is used to evaluate the power of the test and compare it with several other applicable tests that were designed under different contexts. Tolerance intervals within the context of the model are derived, thus generalizing well-known tolerance intervals for ordinary least squares regression. Use of the model and its associated analyses is illustrated with an aerospace application where hundreds of electronic components are continuously monitored by an automated system that flags components that are suspected of unusual degradation patterns.


---

## Body

## 1. Introduction

### 1.1. Background

The model and analyses developed in this paper address a problem encountered when analyzing data from service life tests of aerospace hardware packages. Data for as many as 700 performance metrics per part type are automatically stored during surveillance testing and subsequently input into a software program where, up to this point, an ordinary least squares line based on normal-theory has been routinely fit to the data using time-in-service as the explanatory variable. The software program also outputs tolerance intervals based on the ordinary least squares analysis. Engineers monitoring this process are alerted only to those cases where observations in the scatter plot fall outside the tolerance intervals or the tolerance interval crosses a given limit within some specified future time interval (e.g., 60 months). In cases where the alert suggests an increasing accelerated degradation, a proactive corrective action (e.g., part replacement) may be initiated. In cases where the alert suggests less than expected degradation, a cursory investigation to determine if the part is being utilized properly is initiated. Tolerance intervals are often similarly used for monitoring environmental applications [1].Compared to having engineers individually examine the data from all the combinations of metrics and part types, the automated monitoring and flagging process is quite cost-effective. However, when the variance of the metric increases with time, the tolerance intervals that result from the ordinary least squares analyses fail in the sense that the intervals become too narrow as time increases. Figure1 illustrates this point, where the 111 observations are the performance metric for one type of part. It is evident from the figure that the random errors associated with the regression model are heteroscedastic. Heteroscedasticity is, of course, not unique to our application. It arises routinely in other applications such as economics, behavioral sciences, social sciences, environmental science, and computer vision. The works in [2–6] provide examples within these disciplines, respectively.Figure 1
(Data Set 1) Pointwise 95%-content tolerance intervals with 90% confidence. Dashed lines correspond to an ordinary least squares analysis and solid lines correspond to the estimated weighted least squares analysis.Figure1 also shows two sets of pointwise 95%-content tolerance intervals constructed at the 90% confidence level. The dashed lines represent the tolerance intervals derived from an ordinary least squares analysis, while the solid lines represent the tolerance intervals derived from an estimated weighted least squares analysis that is described in Section 1.2. The bold line in the figure is the estimated weighted least squares regression line. It is evident in Figure 1 that the ordinary least squares tolerance intervals are inadequate in the sense that they are too wide for small time values and too narrow for large time values. On the other hand, the more pronounced curvature associated with the tolerance intervals derived from the estimated weighted least squares analysis more adequately captures the range of the observations at both ends of the time spectrum.


### 1.2. Model for Heteroscedasticity

Examination of aerospace data for many part types and many performance metrics led us to propose the following model for heteroscedasticity. LetX(t) denote the metric level at time t for a particular part, and for motivational purposes, assume that the underlying process is a discrete-time valued process (t = 0,1,…). Assume that the initial value of the process is X(0)=β0+e, where β0 is an unknown constant and e is a normally distributed random variable with zero mean and unknown variance σe2>0. If the degradation process has both a constant deterministic trend, say β1, and also a random stochastic perturbation, it would follow that X(t)=X(t-1)+β1+st, where st is a normally distributed random variable with zero mean and unknown variance σs2≥0. It follows from successive substitutions that X(t)=β0+β1t+e+∑i=1tsi, or equivalently, X(t)=β0+β1t+e+δt, where δt is a normally distributed random variable with zero mean and unknown variance σs2t. The model for X(t) is (drifting) Brownian motion, except for the fact that the variance function is displaced from the origin by σe2. In our application, the process X(t) is measured one time on each unit; so the data are a collection of independent observations Yi=Xi(ti), where Xi(t) is the process associated with the ith unit. Thus, we have the model Yi=β0+β1ti+ei+δi (i=1,…,n) for the observations.Equivalently, the model implies that observations are independent and normally distributed with meansβ0+β1ti and variances σe2+σs2ti. The fact that a variance component, σs2, is responsible for the heteroscedasticity is somewhat unique to our model since the large literature pertaining to heteroscedasticity models usually assumes variances of the form h(σe2+λti), where h(·) is an arbitrary function (e.g., a polynomial or the exponential function); see [2, 7] and references therein. While these types of models have seen many successful applications, the lucid interpretation of the heteroscedastic variance function σe2+σs2t was crucial when eliciting buy-in from our aerospace engineering project stakeholders.


### 1.3. Tolerance Interval for Heteroscedastic Regression

Definingρ=σs2/σe2, the log-likelihood function for (β0,β1,σe2,ρ) based on the observations is(1)l(β0,β1,σe2,ρ)=-n2logσe2-12∑i=1nlog(1+ρti)-12σe2∑i=1n(yi-β0-β1ti)21+ρti.
For a fixed ρ the conditional MLEs, say β̃0(ρ) and β̃1(ρ), of β0 and β1 satisfy the linear system of equations:(2)β̃0(ρ)∑i=1n11+ρti+β̃1(ρ)∑i=1nti1+ρti=∑i=1nyi1+ρti,β̃0(ρ)∑i=1nti1+ρti+β̃1(ρ)∑i=1nti21+ρti=∑i=1nyiti1+ρti,
and the conditional MLE of σe2 is σ̃e2(ρ)=[∑i=1n(yi-β̃0(ρ)-β̃1(ρ)ti)2/(1+ρti)]/n. The conditional MLEs, β̃0(ρ) and β̃1(ρ), are the weighted least squares estimator of the slope and intercept and have an estimated variance-covariance matrix equal to Σ(ρ)=σ̃e2(ρ)[X′Ω-1(ρ)X]-1, where X is the n×2 matrix whose first column is all ones and whose second column is the set of values {ti}i=1n and Ω(ρ) = Diag(1+ρt1,…,1+ρtn).The profile log-likelihood function (e.g., see reference [8]) of ρ is lp(ρ)=-nlogσ̃e2(ρ)/2-∑i=1nlog(1+ρti)/2-n/2, which can be maximized to find the MLE of ρ, say ρ̂. The MLEs of β0, β1 and σe2 are then are β̂0=β̃0(ρ̂), β̂1=β̃1(ρ̂), and σ̂e2=σ̃e2(ρ̂). Unbiasedness of β̂0 and β̂1 follows from general results in [9].A pointwise100γ%-content tolerance interval with confidence level 100(1-α)% for a regression line in the context of homoscedastic normal errors was derived in [10]. The following theorem, which is proved in the Appendix A, extends that result to our context.Theorem 1.
Ifρ is known, a pointwise 100γ%-content tolerance interval with confidence level 100(1-α)% for a normal distribution that has mean β0+β1t and variance σe2(1+ρt) is
(3)β̃0(ρ)+β̃1(ρ)t±r(ρ)σ̃e2(ρ)(1+ρt)nχn-2,1-α2,
where χn-2,1-α2 is the upper 1-α percentile of a chi-square distribution with n-2 degrees of freedom and r(ρ) is the solution to (2π)-1/2∫1/N(ρ)-r(ρ)1/N(ρ)+r(ρ)e-t2/2dt=γ, with N(ρ)=(1+ρt)/[(1,t)(X′Ω-1(ρ)X)-1(1,t)′].Proof.
See AppendixA.The pointwise tolerance intervals shown in Figure1 corresponding to the estimated weighted least squares analysis were drawn from (3) with ρ replaced by its MLE ρ̂, and with α=0.1 and γ=0.95.


### 1.4. Testing for Heteroscedasticity

Of special importance when analyzing the data for our application is the test ofH0 that the data are homoscedastic. If the hypothesis is not rejected, all the pertinent inference can be drawn from an efficient ordinary least squares analyses. On the other hand, rejecting H0 protects a practitioner from the pitfalls of using an inappropriate least squares analysis which include a higher than expected false alert rate. Deriving a test for H0 is particularly interesting from the standpoint that there are alternative tests available in the literature and in subsequent sections in this paper we present a unifying framework for relating alternative tests to each other.


### 1.5. Overview of Remaining Sections

In Section2, we position the proposed heteroscedastic model as a special case of a general mixed linear model. We then propose a test of H0 that is motivated by informally combining the root arguments used when deriving locally most powerful tests and score tests. We show that the test statistic has a distribution that only depends on the underlying variance ratio ρ and thus is amenable to a significance test of H0:ρ=0. A Monte Carlo algorithm for computing the required critical values is outlined. In Section 3, we relate our proposed test to some common tests for heteroscedasticity, namely, the Breusch and Pagan, White, and likelihood ratio tests. In Section 4, we show additional examples of our test for heteroscedasticity and the use of tolerance intervals within the context of aerospace engineering case study. We also report results of a simulation study that examined the power of all the tests that are discussed in this paper. We close the paper with a short summary in Section 5.


## 1.1. Background

The model and analyses developed in this paper address a problem encountered when analyzing data from service life tests of aerospace hardware packages. Data for as many as 700 performance metrics per part type are automatically stored during surveillance testing and subsequently input into a software program where, up to this point, an ordinary least squares line based on normal-theory has been routinely fit to the data using time-in-service as the explanatory variable. The software program also outputs tolerance intervals based on the ordinary least squares analysis. Engineers monitoring this process are alerted only to those cases where observations in the scatter plot fall outside the tolerance intervals or the tolerance interval crosses a given limit within some specified future time interval (e.g., 60 months). In cases where the alert suggests an increasing accelerated degradation, a proactive corrective action (e.g., part replacement) may be initiated. In cases where the alert suggests less than expected degradation, a cursory investigation to determine if the part is being utilized properly is initiated. Tolerance intervals are often similarly used for monitoring environmental applications [1].Compared to having engineers individually examine the data from all the combinations of metrics and part types, the automated monitoring and flagging process is quite cost-effective. However, when the variance of the metric increases with time, the tolerance intervals that result from the ordinary least squares analyses fail in the sense that the intervals become too narrow as time increases. Figure1 illustrates this point, where the 111 observations are the performance metric for one type of part. It is evident from the figure that the random errors associated with the regression model are heteroscedastic. Heteroscedasticity is, of course, not unique to our application. It arises routinely in other applications such as economics, behavioral sciences, social sciences, environmental science, and computer vision. The works in [2–6] provide examples within these disciplines, respectively.Figure 1
(Data Set 1) Pointwise 95%-content tolerance intervals with 90% confidence. Dashed lines correspond to an ordinary least squares analysis and solid lines correspond to the estimated weighted least squares analysis.Figure1 also shows two sets of pointwise 95%-content tolerance intervals constructed at the 90% confidence level. The dashed lines represent the tolerance intervals derived from an ordinary least squares analysis, while the solid lines represent the tolerance intervals derived from an estimated weighted least squares analysis that is described in Section 1.2. The bold line in the figure is the estimated weighted least squares regression line. It is evident in Figure 1 that the ordinary least squares tolerance intervals are inadequate in the sense that they are too wide for small time values and too narrow for large time values. On the other hand, the more pronounced curvature associated with the tolerance intervals derived from the estimated weighted least squares analysis more adequately captures the range of the observations at both ends of the time spectrum.


## 1.2. Model for Heteroscedasticity

Examination of aerospace data for many part types and many performance metrics led us to propose the following model for heteroscedasticity. LetX(t) denote the metric level at time t for a particular part, and for motivational purposes, assume that the underlying process is a discrete-time valued process (t = 0,1,…). Assume that the initial value of the process is X(0)=β0+e, where β0 is an unknown constant and e is a normally distributed random variable with zero mean and unknown variance σe2>0. If the degradation process has both a constant deterministic trend, say β1, and also a random stochastic perturbation, it would follow that X(t)=X(t-1)+β1+st, where st is a normally distributed random variable with zero mean and unknown variance σs2≥0. It follows from successive substitutions that X(t)=β0+β1t+e+∑i=1tsi, or equivalently, X(t)=β0+β1t+e+δt, where δt is a normally distributed random variable with zero mean and unknown variance σs2t. The model for X(t) is (drifting) Brownian motion, except for the fact that the variance function is displaced from the origin by σe2. In our application, the process X(t) is measured one time on each unit; so the data are a collection of independent observations Yi=Xi(ti), where Xi(t) is the process associated with the ith unit. Thus, we have the model Yi=β0+β1ti+ei+δi (i=1,…,n) for the observations.Equivalently, the model implies that observations are independent and normally distributed with meansβ0+β1ti and variances σe2+σs2ti. The fact that a variance component, σs2, is responsible for the heteroscedasticity is somewhat unique to our model since the large literature pertaining to heteroscedasticity models usually assumes variances of the form h(σe2+λti), where h(·) is an arbitrary function (e.g., a polynomial or the exponential function); see [2, 7] and references therein. While these types of models have seen many successful applications, the lucid interpretation of the heteroscedastic variance function σe2+σs2t was crucial when eliciting buy-in from our aerospace engineering project stakeholders.


## 1.3. Tolerance Interval for Heteroscedastic Regression

Definingρ=σs2/σe2, the log-likelihood function for (β0,β1,σe2,ρ) based on the observations is(1)l(β0,β1,σe2,ρ)=-n2logσe2-12∑i=1nlog(1+ρti)-12σe2∑i=1n(yi-β0-β1ti)21+ρti.
For a fixed ρ the conditional MLEs, say β̃0(ρ) and β̃1(ρ), of β0 and β1 satisfy the linear system of equations:(2)β̃0(ρ)∑i=1n11+ρti+β̃1(ρ)∑i=1nti1+ρti=∑i=1nyi1+ρti,β̃0(ρ)∑i=1nti1+ρti+β̃1(ρ)∑i=1nti21+ρti=∑i=1nyiti1+ρti,
and the conditional MLE of σe2 is σ̃e2(ρ)=[∑i=1n(yi-β̃0(ρ)-β̃1(ρ)ti)2/(1+ρti)]/n. The conditional MLEs, β̃0(ρ) and β̃1(ρ), are the weighted least squares estimator of the slope and intercept and have an estimated variance-covariance matrix equal to Σ(ρ)=σ̃e2(ρ)[X′Ω-1(ρ)X]-1, where X is the n×2 matrix whose first column is all ones and whose second column is the set of values {ti}i=1n and Ω(ρ) = Diag(1+ρt1,…,1+ρtn).The profile log-likelihood function (e.g., see reference [8]) of ρ is lp(ρ)=-nlogσ̃e2(ρ)/2-∑i=1nlog(1+ρti)/2-n/2, which can be maximized to find the MLE of ρ, say ρ̂. The MLEs of β0, β1 and σe2 are then are β̂0=β̃0(ρ̂), β̂1=β̃1(ρ̂), and σ̂e2=σ̃e2(ρ̂). Unbiasedness of β̂0 and β̂1 follows from general results in [9].A pointwise100γ%-content tolerance interval with confidence level 100(1-α)% for a regression line in the context of homoscedastic normal errors was derived in [10]. The following theorem, which is proved in the Appendix A, extends that result to our context.Theorem 1.
Ifρ is known, a pointwise 100γ%-content tolerance interval with confidence level 100(1-α)% for a normal distribution that has mean β0+β1t and variance σe2(1+ρt) is
(3)β̃0(ρ)+β̃1(ρ)t±r(ρ)σ̃e2(ρ)(1+ρt)nχn-2,1-α2,
where χn-2,1-α2 is the upper 1-α percentile of a chi-square distribution with n-2 degrees of freedom and r(ρ) is the solution to (2π)-1/2∫1/N(ρ)-r(ρ)1/N(ρ)+r(ρ)e-t2/2dt=γ, with N(ρ)=(1+ρt)/[(1,t)(X′Ω-1(ρ)X)-1(1,t)′].Proof.
See AppendixA.The pointwise tolerance intervals shown in Figure1 corresponding to the estimated weighted least squares analysis were drawn from (3) with ρ replaced by its MLE ρ̂, and with α=0.1 and γ=0.95.


## 1.4. Testing for Heteroscedasticity

Of special importance when analyzing the data for our application is the test ofH0 that the data are homoscedastic. If the hypothesis is not rejected, all the pertinent inference can be drawn from an efficient ordinary least squares analyses. On the other hand, rejecting H0 protects a practitioner from the pitfalls of using an inappropriate least squares analysis which include a higher than expected false alert rate. Deriving a test for H0 is particularly interesting from the standpoint that there are alternative tests available in the literature and in subsequent sections in this paper we present a unifying framework for relating alternative tests to each other.


## 1.5. Overview of Remaining Sections

In Section2, we position the proposed heteroscedastic model as a special case of a general mixed linear model. We then propose a test of H0 that is motivated by informally combining the root arguments used when deriving locally most powerful tests and score tests. We show that the test statistic has a distribution that only depends on the underlying variance ratio ρ and thus is amenable to a significance test of H0:ρ=0. A Monte Carlo algorithm for computing the required critical values is outlined. In Section 3, we relate our proposed test to some common tests for heteroscedasticity, namely, the Breusch and Pagan, White, and likelihood ratio tests. In Section 4, we show additional examples of our test for heteroscedasticity and the use of tolerance intervals within the context of aerospace engineering case study. We also report results of a simulation study that examined the power of all the tests that are discussed in this paper. We close the paper with a short summary in Section 5.


## 2. Proposed Test for Heteroscedasticity

The heteroscedastic regression model developed in Section1 is a special case of a more general mixed linear model defined by(4)Y=Xβ+Zs+e,
where X is an n×p matrix (we assume without loss of generality that X has a full column rank) of known constants, β is p×1 vector of fixed-effects, Z is n×m matrix of known constants, s is a m×1 vector of unobservable random-effects that follow an MVNm(0,σs2Im) distribution, and e is an n×1 vector of unobservable random error terms, independent of s, and having a MVNn(0,σe2In) distribution. The special case considered in Section 1 is simply where p=2 and m=n with Z a diagonal matrix with the values equal to {ti}i=1n.We begin the derivation of our proposed test for heteroscedasticity by writing the log-likelihood of the general mixed linear model, based on(β,σe2,ρ), as follows:(5)l(β,σe2,ρ)∝-n(logσe2)2-log|I+ρZZ′|2-(y-Xβ)′(I+ρZZ′)-1(y-Xβ)2σe2.
Next, there exists an orthogonal matrix Q such that Q′ZZ′Q=B, where B is a diagonal matrix whose elements {λi}i=1n are the eigenvalues of ZZ′. It follows that ZZ′=QBQ′, and hence Q′(I+ρZZ′)Q=I+ρB. Equivalently, Q′(I+ρZZ′)-1Q=(I+ρB)-1. Combining this result with the definitions v=Q′y and μ=Q′Xβ, the log-likelihood can be rewritten as(6)l(β,σe2,ρ)∝-n2logσe2-12log|I+ρB|-12σe2(Q′y-Q′Xβ)′(I+ρB)-1(Q′y-Q′Xβ)=-n2logσe2-12∑i=1nlog(1+ρλi)-12σe2∑i=1n(vi-μi)21+ρλi.
If β and σe2 were known, the locally most powerful test statistic would be(7)∂logL(β,σe2,ρ)∂ρ|ρ=0∝∑i=1n(vi-μi)2λiσe2=(Q′y-Q′Xβ)′B(Q′y-Q′Xβ)σe2.
To obtain a useable test statistic we follow the procedure used when deriving the score test (e.g., see [2]) by replacing the unknown β and σe2 by their conditional MLEs given H0, which are β̃0(0)=(X′X)-1X′y and σ̃e2(0)=r′r/n, where r=(I-PX)y and PX=X(X′X)-1X′. The resulting statistic is(8)R=(Q′y-Q′PXy)′B(Q′y-Q′PXy)r′r=r′ZZ′rr′r.The distribution ofr is MVNn{0,σe2[I-PX+ρ(I-PX)ZZ′(I-PX)]} and by the nature of R it is clear that its distribution only depends on ρ. A simple approach for finding critical values to test H0 is to simulate r vectors from an MVNn(0,I-PX), compute the corresponding values of R from (8), and estimate the critical value from the resulting empirical distribution (edf) function of R values. Denoting the upper-α percentile of the edf by Ra, the R test rejects H0 if R>Rα. We emphasize the computational simplicity of R, and the ease in which critical values can be obtained. It is also possible to reexpress the null distribution of (8) as a ratio of quadratic forms in independent standard normal random variables and then the saddlepoint approximation described by [11] can be employed to approximate the P-value associated with the test.A locally most powerful invariant (LMPI) test ofH0 was derived in [12–14]. The derivation of this test is sketched in Appendix B. It can be seen there that the complexity of computing LMPI test is nontrivial. The key result in Appendix B is the demonstration that the LMPI test and the R test are equivalent under the general mixed linear model. The practical significance of this finding is that the computation of the LMPI test, as originally proposed, involves a complicated eigenvalue and eigenvector analysis, whereas computing the equivalent form R test only requires fitting the model with ordinary least squares tools. To characterize this computational contrast another way, the computations associated with the equivalent form R test can be done with software as simple as Excel, whereas the computations associated with the originally proposed form LMPI test will require a rich matrix processing software package.From a conceptual perspective, the equivalence between the LMPI test and theR test is interesting in the sense that the original complicated derivation of the LMPI test can be replaced by the simpler derivation of the R test. In addition, the demonstrated equivalence also provides further motivation for a procedure introduced in [15] where a test statistic is derived by replacing nuisance parameters in the locally most powerful test statistic by their conditional MLEs. Results in [15] show this procedure produces an asymptotically uniformly most powerful test. Our equivalence finding shows that when we used this approach to derive the R test, in a finite sample context, we arrive at a locally most powerful invariant test.


## 3. Other Tests for Heteroscedasticity

The Breusch and Pagan test for homogeneity is a partial score test derived for a general setting where the observations are independent and normally distributed with meansxi′β and variances σi2=h(zi′ξ), where xi′ is a p×1 vector of known covariates, β is a p×1 vector of unknown parameters, zi′ is a q×1 vector of known covariates (whose components may or may not overlap with the components of xi′), ξ is a q×1 vector of unknown parameters, and h(·) is an arbitrary function that is only assumed to possess a first derivative. By setting the first component of zi′ to unity, a test of H0*:ξ2=⋯=ξq=0 is a test of heteroscedasticity. The Breusch and Pagan test [16] is the partial score test of H0*. For application to our problem, we have xi=zi=(1,ti)′, ξ=(σe2,σs2)′ and h(zi′ξ)=zi′ξ and the test statistic simplifies to(9)BP=(r′Arr′r)2,
where A is an n×n diagonal matrix whose ith diagonal element is equal to n(ti-t̅)/2∑i=1n(ti-t̅)2.The White test [17] is positioned to be robust to the normality assumption. For our application, the statistic can be obtained by first computing the vector of squared residuals, say u where ui=ri2, and then computing the test statistic as(10)W=nu′(PV-Jn/n)uu′(In-Jn/n)u,
where Jn is an n×n matrix of ones, PV=V(V′V)-1V′, and V is a matrix whose first column is all ones, whose second column is the values {ti}i=1n and whose third column is the values {ti2}i=1n. Note that the ratio term in the definition of W is simply the R-square value of the regression of u on the V matrix. The distributions of both BP and W depend on ρ and the same Monte Carlo method for approximating the null distribution of R can be used to approximate the null distributions of BP and W.The log-likelihood function for the general mixed linear model was given in (5) and computation of the MLE of (β,σe2,ρ) and the conditional MLE of (β,σe2) was discussed in Section 2. It follows that the likelihood ratio test statistic is(11)LRT=-2{l(β̃,σ̃e2,0)-l(β̂,σ̂e2,ρ̂)}.
The LRT requires relatively significant computations. For example, computations of the full model MLEs needed by the LRT are not readily available unless software such as PROC MIXED in SAS is available. In standard situations the asymptotic null distribution of LRT is chi-square with degrees of freedom equal to the difference in the dimension of the unconstrained and constrained parameter spaces. However, our situation is not standard since the reduced parameter space lies on the boundary of the unconstrained parameter space. Results in [18] show that the asymptotic distribution of LRT in (11) is a 50: 50 mixture of zero and a chi-square distribution with one degree of freedom. A size-α LRT therefore rejects H0:ρ=0 if LRT>χ1,2α2, where χ1,2α2 denotes the upper-2α percentile of a chi-square distribution with one degree of freedom.


## 4. Application to Aerospace Case Study

### 4.1. Illustrative Examples

Figure1 was used in Section 1 as illustrative of a data set that exhibits heteroscedasticity. Figure 2 shows a second illustrative data set (109 observations) which does not exhibit evidence of heteroscedasticity. As in Figure 1, two sets of pointwise 95%-content tolerance intervals constructed at the 90% confidence level are shown corresponding to an ordinary least squares analysis and the estimated weighted least squares analysis described in Section 1.2. In contrast to Figure 1, the two sets of tolerance intervals are nearly identical.Figure 2
(Data Set 2) Pointwise 95%-content tolerance intervals with 90% confidence. Dashed lines correspond to an ordinary least squares analysis and solid lines correspond to the estimated weighted least squares analysis.Table1 shows the computed test statistics R, BP, W, and LRT for the two data sets shown in Figures 1 and 2, along with their P-values for testing H0:ρ=0. The tests all agree that heteroscedasticity appears to exist in Data Set 1 but not in Data Set 2. We note that of all the tests considered, R is by far the simplest to compute with BP and W only slightly more complicated. As described in the introduction, our application context has hundreds of metrics that need to be monitored with this type of analysis methodology. Figures 3, 4, 5, and 6 provide four additional illustrative examples, including the computed R statistic and associated P-value for H0:ρ=0. There is evidence of heteroscedasticity in all four of these data sets.Table 1
Tests for Heteroscedasticity in Figures1 and 2.
TestData SetData Set 1Data Set 2StatisticP-valueStatisticP-valueR8.92<.0014.83.062BP52.1<.0012.05.15W26.6<.0013.29.17LRT34.8<.0011.59.10Figure 3
(Data Set 3) Pointwise 95% tolerance intervals with 90% confidence.Figure 4
(Data Set 4) Pointwise 95% tolerance intervals with 90% confidence.Figure 5
(Data Set 5) Pointwise 95% tolerance intervals with 90% confidence.Figure 6
(Data Set 6) Pointwise 95% tolerance intervals with 90% confidence.


### 4.2. Power Comparison of Alternative Tests

It is easy to show that the distributions of theR, BP, W, and LRT test statistics depend only on ρ. It follows that the power functions (i.e., the probability of rejecting H0:ρ=0 as function of ρ) of the tests can be evaluated using simulated data sets from the mixed model in (3) using the parameters (β0,β1,σe2,ρ)=(0,0,1,ρ). We implemented our simulation study using the S-Plus package. For each simulated data set, each of the test statistics can be evaluated and compared to their respective critical values corresponding to nominal size α significance tests. The power of the tests is estimated by the fraction of data sets for which they reject H0.Column 2 of Table2 shows the estimated power of the 10% nominal size R test based on this procedure for ρ∈[0,0.25], using 5000 simulated data sets and the same set of {ti}i=1n values associated with Data Set 1, which is shown in Figure 1. Columns 3–5 of Table 2 give the ratio of the power for the LRT, BP, and W tests, relative to the R test. The results show that the R test and the LRT have nearly identical power, though the R test has slightly better power for small ρ. The R test has uniformly better power than BP and W tests. A somewhat surprising observation is the pronounced dominance of R with respect to BP and W, which are perhaps the two most widely known tests in the economics literature. Table 3 shows similar power comparisons using the {ti}i=1n associated with Data Set 2, that is shown in Figure 2, and the conclusions are consistent with those drawn from Table 2.Table 2
Power comparison of the tests using{ti}i=1n values from Data Set 1.
ρPower ofR TestRatio of Power Relative toR TestLRTBPW0.000.1000.940.910.920.010.1550.970.670.670.020.2060.970.620.640.030.2781.000.630.650.040.3440.960.660.600.050.3920.990.680.620.060.4431.020.680.670.070.4921.010.720.670.080.5411.010.720.680.090.5841.010.770.710.100.6241.020.770.700.110.6631.020.770.730.120.7011.010.800.730.130.7191.020.800.730.140.7411.030.820.760.150.7631.010.840.770.160.7831.020.840.780.170.8091.010.850.780.180.8301.010.860.800.190.8301.020.880.820.200.8431.020.870.840.210.8611.010.880.830.220.8691.010.880.840.230.8751.020.900.840.240.8841.020.900.850.250.9001.020.900.850.300.9231.010.930.890.350.9431.000.940.920.400.9491.010.960.93Table 3
Power comparison of the tests using{ti}i=1n values from Data Set 2.
ρPower ofR TestRatio of Power Relative toR TestLRTBPW00.1001.010.960.890.0010.1070.950.910.840.0020.1041.010.900.980.0030.1260.920.780.750.0040.1270.960.720.810.0050.1280.920.750.680.0060.1341.000.700.810.0070.1460.840.650.650.0080.1361.030.780.760.0090.1371.060.780.770.0100.1480.900.730.680.0110.1521.010.720.780.0120.1740.860.650.680.0130.1710.970.680.730.0140.1681.030.680.740.0150.1830.920.630.700.0160.1890.930.690.670.0170.2000.920.630.650.0180.2040.920.650.620.0190.2050.990.670.650.020.2060.990.640.640.030.2670.980.660.610.040.3161.100.660.640.050.3840.850.710.630.060.4431.010.710.630.070.4961.020.740.640.080.5441.020.740.680.090.5861.030.760.700.100.6391.000.760.680.110.6761.010.800.710.120.7001.010.820.720.130.7351.010.830.770.140.7561.020.840.770.150.7831.010.840.790.160.7991.010.870.810.170.8281.000.850.800.180.8301.030.870.820.190.8531.010.890.830.200.8681.010.890.840.210.8781.010.900.850.220.8871.010.910.860.230.8971.010.910.860.240.9031.020.920.890.250.9071.030.930.890.300.9351.020.950.920.350.9591.010.960.940.400.9731.010.960.95


## 4.1. Illustrative Examples

Figure1 was used in Section 1 as illustrative of a data set that exhibits heteroscedasticity. Figure 2 shows a second illustrative data set (109 observations) which does not exhibit evidence of heteroscedasticity. As in Figure 1, two sets of pointwise 95%-content tolerance intervals constructed at the 90% confidence level are shown corresponding to an ordinary least squares analysis and the estimated weighted least squares analysis described in Section 1.2. In contrast to Figure 1, the two sets of tolerance intervals are nearly identical.Figure 2
(Data Set 2) Pointwise 95%-content tolerance intervals with 90% confidence. Dashed lines correspond to an ordinary least squares analysis and solid lines correspond to the estimated weighted least squares analysis.Table1 shows the computed test statistics R, BP, W, and LRT for the two data sets shown in Figures 1 and 2, along with their P-values for testing H0:ρ=0. The tests all agree that heteroscedasticity appears to exist in Data Set 1 but not in Data Set 2. We note that of all the tests considered, R is by far the simplest to compute with BP and W only slightly more complicated. As described in the introduction, our application context has hundreds of metrics that need to be monitored with this type of analysis methodology. Figures 3, 4, 5, and 6 provide four additional illustrative examples, including the computed R statistic and associated P-value for H0:ρ=0. There is evidence of heteroscedasticity in all four of these data sets.Table 1
Tests for Heteroscedasticity in Figures1 and 2.
TestData SetData Set 1Data Set 2StatisticP-valueStatisticP-valueR8.92<.0014.83.062BP52.1<.0012.05.15W26.6<.0013.29.17LRT34.8<.0011.59.10Figure 3
(Data Set 3) Pointwise 95% tolerance intervals with 90% confidence.Figure 4
(Data Set 4) Pointwise 95% tolerance intervals with 90% confidence.Figure 5
(Data Set 5) Pointwise 95% tolerance intervals with 90% confidence.Figure 6
(Data Set 6) Pointwise 95% tolerance intervals with 90% confidence.


## 4.2. Power Comparison of Alternative Tests

It is easy to show that the distributions of theR, BP, W, and LRT test statistics depend only on ρ. It follows that the power functions (i.e., the probability of rejecting H0:ρ=0 as function of ρ) of the tests can be evaluated using simulated data sets from the mixed model in (3) using the parameters (β0,β1,σe2,ρ)=(0,0,1,ρ). We implemented our simulation study using the S-Plus package. For each simulated data set, each of the test statistics can be evaluated and compared to their respective critical values corresponding to nominal size α significance tests. The power of the tests is estimated by the fraction of data sets for which they reject H0.Column 2 of Table2 shows the estimated power of the 10% nominal size R test based on this procedure for ρ∈[0,0.25], using 5000 simulated data sets and the same set of {ti}i=1n values associated with Data Set 1, which is shown in Figure 1. Columns 3–5 of Table 2 give the ratio of the power for the LRT, BP, and W tests, relative to the R test. The results show that the R test and the LRT have nearly identical power, though the R test has slightly better power for small ρ. The R test has uniformly better power than BP and W tests. A somewhat surprising observation is the pronounced dominance of R with respect to BP and W, which are perhaps the two most widely known tests in the economics literature. Table 3 shows similar power comparisons using the {ti}i=1n associated with Data Set 2, that is shown in Figure 2, and the conclusions are consistent with those drawn from Table 2.Table 2
Power comparison of the tests using{ti}i=1n values from Data Set 1.
ρPower ofR TestRatio of Power Relative toR TestLRTBPW0.000.1000.940.910.920.010.1550.970.670.670.020.2060.970.620.640.030.2781.000.630.650.040.3440.960.660.600.050.3920.990.680.620.060.4431.020.680.670.070.4921.010.720.670.080.5411.010.720.680.090.5841.010.770.710.100.6241.020.770.700.110.6631.020.770.730.120.7011.010.800.730.130.7191.020.800.730.140.7411.030.820.760.150.7631.010.840.770.160.7831.020.840.780.170.8091.010.850.780.180.8301.010.860.800.190.8301.020.880.820.200.8431.020.870.840.210.8611.010.880.830.220.8691.010.880.840.230.8751.020.900.840.240.8841.020.900.850.250.9001.020.900.850.300.9231.010.930.890.350.9431.000.940.920.400.9491.010.960.93Table 3
Power comparison of the tests using{ti}i=1n values from Data Set 2.
ρPower ofR TestRatio of Power Relative toR TestLRTBPW00.1001.010.960.890.0010.1070.950.910.840.0020.1041.010.900.980.0030.1260.920.780.750.0040.1270.960.720.810.0050.1280.920.750.680.0060.1341.000.700.810.0070.1460.840.650.650.0080.1361.030.780.760.0090.1371.060.780.770.0100.1480.900.730.680.0110.1521.010.720.780.0120.1740.860.650.680.0130.1710.970.680.730.0140.1681.030.680.740.0150.1830.920.630.700.0160.1890.930.690.670.0170.2000.920.630.650.0180.2040.920.650.620.0190.2050.990.670.650.020.2060.990.640.640.030.2670.980.660.610.040.3161.100.660.640.050.3840.850.710.630.060.4431.010.710.630.070.4961.020.740.640.080.5441.020.740.680.090.5861.030.760.700.100.6391.000.760.680.110.6761.010.800.710.120.7001.010.820.720.130.7351.010.830.770.140.7561.020.840.770.150.7831.010.840.790.160.7991.010.870.810.170.8281.000.850.800.180.8301.030.870.820.190.8531.010.890.830.200.8681.010.890.840.210.8781.010.900.850.220.8871.010.910.860.230.8971.010.910.860.240.9031.020.920.890.250.9071.030.930.890.300.9351.020.950.920.350.9591.010.960.940.400.9731.010.960.95


## 5. Summary

Our aerospace case study, which pertains to the time evolution of performance metrics for electronic equipment, motivated the derivation of a model for heteroscedastic regression errors. A test for homoscedasticity was proposed and compared with various other tests that are prevalent in the literature. The proposedR test was shown to be equivalent to an LMPI test associated with a general mixed linear model of which this paper’s heteroscedastic regression model can be regarded as a special case. Aside from what we judge to be a simpler intuitive motivation associated with the R test formulation, it has the practical benefit of being much simpler to compute than the formula that results from the classical LMPI derivation. The power of the R test was compared to some popular alternative tests, and only the considerably more difficult to compute likelihood ratio test has power that is comparable to the R test. We extended the classical application of tolerance intervals for regression to our heteroscedastic case and illustrated their use with six data sets from our aerospace case study.



---

*Source: 126283-2010-02-28.xml*