Likelihood Ratio Test and Non-parametric Test for Load Sharing

In present article, we propose a likelihood ratio test and a non-parametric test for testing the load sharing effect observed in the two component parallel load sharing system. We have modeled the load sharing phenomenon observed in such system by the exponentiated conditional distribution function based load sharing model proposed by Sutar and Naik-Nimbalkar (2016). We have taken component baseline distribution as Weibull distribution and linear failure rate distribution. The simulation study to see the performance of proposed test procedures is reported. We analyze two data sets for illustrative purpose.


Introduction
The systems in which the failure of a component increases or decreases the load on other surviving components, thereby decreasing or increasing their chances of survival are known as load sharing systems. Such systems are also referred to as dynamic reliability systems and have several practical applications. Examples of load sharing systems include fibrous composite materials, power plants, automobiles, two jet engines of an airplane among others (Liu (1998)). These load sharing models have many applications in various fields engineering (Singh and Gupta (2012)). some of the references cited there in. Sutar and Naik-Nimbalkar (2016) have considered a k-out-of-m system and proposed load sharing model through exponentiated conditional distributions of the order statistics. It leads to a new family of multivariate distributions for the ordered r.v.'s, which is different from the one introduced by Kamps (1995) and Cramer and Kamps (1996), Cramer and Kamps (2001) for sequential order statistics. In the present article, we consider a model for k-out-of-m load sharing systems proposed by Sutar and Naik-Nimbalkar (2016). We consider a two component parallel load sharing system and propose the two tests namely the likelihood ratio test (LRT) and the non-parametric test for testing the load sharing effect observed in such system. The performance of the proposed test procedures is also evaluated with extensive simulation study. The rest of this article is organized as follows.
In Section 2, we restate the load sharing model for a two component parallel system, in a general setup and also when the initial lifetimes of the components having Weibull and linear failure rate distributions. For the general setup, we propose a LRT and non-parametric tests to test the hypothesis that, the component failures occurs independently against the alternative that they show the load sharing phenomenon in Section 3. In Section 4, simulation study showing the performance of the test procedures is reported. We report the illustrations in Section 5 and conclusions are given in the last section.

The model and the associated parameter estimation procedure 2.1. General setup
We consider a two component parallel system with U 1 , U 2 as the component lifetimes, being independent and identically distributed (i.i.d.) with common probability density function (p.d.f.) f (·; λ) where λ may be scalar or vector valued parameter. Let F (·; λ) andF (·; λ) denote corresponding cumulative distribution function (c.d.f.) and survival function (s.f.). Let X = min(U 1 , U 2 ) denote the first failure time and Y = max(U 1 , U 2 ) denote the second failure time. The density of X is k X (x, λ) = 2f (x; λ)F (x; λ), x > 0.
According to load sharing model proposed by Sutar and Naik-Nimbalkar (2016), the conditional density of Y given X = x is and the joint density of (X, Y ) is We note that, β = 1 implies independent setup whereas β = 1 implies dependent or load sharing setup. Let {(x i , y i ), i = 1, 2, · · · , n} be a random sample of size n from the joint density given in (2). We use the two step procedure proposed by Sutar and Naik-Nimbalkar (2016) for estimating the unknown parameters. This procedure is described as follow.
Step (i): Estimate λ based on observations on X, the first failure time.
That is we can estimate λ by maximizing the following log-likelihood, log L X (λ),based on X, Step (ii): Estimate β by using conditional distribution of Y given X = x and by plugging in the estimates of λ asλ (say), obtained in step (i). Thus, the estimateβ of β is given bŷ In the next subsection, we discuss the same setup when baseline distributions are Weibull and LFR.

Weibull distribution
We consider a two component parallel load sharing system, with the component lifetimes being i.i.d. Weibull r.v.'s with shape parameter λ and scale parameter θ, with p.d.f.
Then, the first failure time, X is a Weibull random variable with shape parameter λ and scale parameter 2θ. The conditional density of Y given X = x as The joint density of (X,Y) can be written as Estimation procedure We adopt two step estimation procedure proposed by Sutar and Naik-Nimbalkar (2016) and the two step estimation procedure for estimating (λ, θ, β) is as follow.
Step (i): Estimation of λ and θ The log-likelihood based on X as log L X = n log 2 + n log θ + n log λ + (λ − 1) The maximum likelihood estimates (λ,θ) of (λ, θ) by solving the following log-likelihood equations by iterative procedures such as Newton-Raphson method, Step (ii): Estimation of β Once we get estimates of λ and θ, the estimate of β follows from the expression (4), thusβ

LFR case
We consider a two component parallel load sharing system, with the component lifetimes being i.i.d. LFR r.v.'s with parameters λ 1 , λ 2 , with p.d.f.
Step (ii): Estimation of β We estimate β from the expression (4) aŝ In the next section, we propose a test procedures for testing the load sharing effect.

The proposed test procedures
We consider a two component parallel load sharing system and our interest is to test whether the lifetimes of the system components are independent against the alternative that there exits a load sharing phenomenon. We propose a LRT and non-parametric test for testing the load sharing effect. We set the null hypothesis as H 0 : β = 1, i.e. the first failure does not affect the lifetime of the second component and the alternative hypothesis as H 1 : β = 1, i.e. the first failure affects (stochastically changes) the lifetime of the second component. The data consist of n i.i.d. pairs of ordered component lifetimes and is denote by, (x, y) = {(x i , y i ) : x i ≤ y i ; i = 1, 2, ..., n}. The joint density of (X, Y ) in terms of conditional of Y given X = x and marginal density of X can be written as Then the complete likelihood, L (X,Y ) , based on data (x, y) is given by , is the likelihood based on marginal density of X.
In the next subsections, we propose two test procedures for testing the load sharing effect.

Likelihood ratio test
The LRT statistic for testing H 0 against H 1 is We modify the usual LRT statistic by estimating β and λ using the two step procedure discussed in the previous section. Let these estimators be denoted byβ andλ. We also note that, sup x)} as the marginal density of X remains the same in both H 0 and H 1 setups. Thus, the test statistic can be written as Lemma 1.λ is consistent estimator of λ 0 (true parameter). That is, as n → ∞, Proof. The consistency ofλ follows from the Cramer-Huzurbazar theorem (Kale (2005)). As the all regularity conditions of the Cramer-Huzurbazar theorem is satisfied by the density k(x, λ) of, the first failure, X. We claim that k(x) = 2f (x; λ)F (x; λ) belongs to the Cramer class.
Asλ is consistent for λ and λ * = bλ 0 + (1 − b)λ, b ∈ (0, 1), we have λ * p → λ 0 . Also we know that if m(·) is continuous function then, by invariance property of consistent estimator, we have m(λ * ) p → m(λ 0 ). Therefore, we can claim that, 1 is bounded in probability. Thus, by using above arguments and from Lehmma 1, we have (β − β 0 ) where Proof. Asβ p → β 0 and β 0 > 0 , we have the result Lemma 4. Asymptotic normality ofβ. That is Proof. By using Taylor series expansion, we expand By the central limit theorem (CLT), we have By Slutsky's theorem on convergence of random variables, we have the following result Theorem 1. Under the Cramer regularity conditions (Cramer (1946)), the asymptotic null distribution of the test statistic, −2 log µ(X, Y) is χ 2 distribution with 1 degree of freedom.
We assume that the joint density k(x, y) of (X, Y ) satisfies the usual regularity conditions. For shake of simplicity, we take λ = λ i.e. scalar parameter and note that the same steps can be extended to λ. The conditional density of Y given X = x is We estimate λ from the marginal distribution of X and β from conditional distribution of Y given X = x by plugging the estimateλ of λ.
From Lemma 2, Lemma 3 and the fact that 1 (β * ) 3 < ∞, we can claim that R n p → 0. Thus, by Lemma 1 to 4, we have the result, that is, as n → ∞, −2 log µ(X, Y ) has the χ 2 distribution with 1 degree of freedom.

Non-parametric test
In this sub-section, we propose the non-parametric test, for testing H 0 : β = 1 against H 1 : β = 1. The non-parametric estimate β of β is obtained by replacingF (·) by its empirical survival function,F n , in expression (4) and is given by We note that, we do not have observations from F (·), baseline distribution, but we do have observation from the minimum (first failure), X. Suppose the c.d.f. of the minimum is K(·) = 1 − (1 − F (·)) 2 , then by this assumption, we can write F (·) = 1 − (1 − K(·)). Let K n denote the empirical distribution function based on observations from K(·), then F n (·) can be estimated as 1 − (1 − K n (·)).
→F (x) and hence, we can show that, where, β is as given in (4).
To obtain asymptotic distribution of β, we define, Z i = − log 1 −F (y i ) F (x i ) , ∀ i = 1, 2, ..., n. Since, (X i , Y i ) being i.i.d. pairs, we can claim that, Z i are also i.i.d. random variable. By CLT, we have,Z By definition, After putting 1 −F (y i ) F (x i ) = u, we get, Using similar argument, we can show that, E(Z 2 1 ) = 2 β 2 . Hence, the variance of Z 1 , V (Z 1 ) is = 1 β 2 . Thus, we claim thatZ By Delta method, we have, Thus, by law of convergence, we have the result.

Simulation study
In the present section, we carry out a simulation study to examine the power and the level attained by the LRT discussed in the Sub-section 3.1. In each case 10000 samples were generated for the different combinations of the sample sizes (n) and the parameter values from the respective densities. There is evidence against H 0 , if the value of the test statistic is ≥ 3.84, which is the critical value at 0.05 of the χ 2 distribution with 1 degree of freedom. In order to estimate power of LRT in case of Weibull, we generate 10000 samples of sample sizes, n = 10, 30 and 50 with parameters (λ, θ) as (1,1), (1,2), (2,1), (0.5,1), (0.5, 2) and β as 0.5, 0.8, 1.5, 1.8 from the density given in (7). To see the level attained by the test, we generated 10000 samples of sample sizes, n = 30, 50, 100 with parameters (λ, θ) as (1,1), (1,2), (2,1), (0.5, 1), (0.5, 2) and β = 1. We report the level attained in Table 1 and the estimated power of LRT for Weibull in Table 2. From Table 2, it is observed that the power improves as the sample size increases and as the parameter β moves away from 1. But we need a large sample size to attain 0.05 level. Similarly, in order to estimate power of LRT in case of LFR, we generate 10000 samples of sample sizes, n = 10, 30 and 50 with parameters (λ 1 , λ 2 ) as (1,1), (1,2), (2,1) and β as 0.5, 0.8, 1.2, 1.5, 1.8 from the density given in (10). To examine the level attained by the test, we generate 10000 samples of sample sizes, = 30, 50, 200, 300, 500 with parameters (λ 1 , λ 2 ) as (1,1), (1,2), (2,1) and β = 1. We report the level attained in Table 3 and the estimated power of LRT for LFR in Table 4. From Table 4, it is observed that the power improves as the sample size increases and as the parameter β moves away from 1. To achieve 0.05, the desired level of significance, we need a large sample size.
The simulation study of the non-parametric test, discussed in the Sub-section 3.2, is also presented. To see the estimated power and level attained, we generate 10000 samples for different combinations of sample sizes and parameter values from the respective densities. Under H 0 : β = 1, the distribution of test statistic √ n( β − 1) for testing H 0 : β = 1 against H 1 : β = 1 is N (0, 1) and we use standard normal cutoff for the proposed non-parametric test procedure. Thus, we reject H 0 , if the calculated value of test statistic ≥ 1.96, the critical value at 5% of standard normal variate. We have observations from K(·) = 1 − (1 − F (·)) m , the distribution function of first failure. Thus, we replace F (·) by 1 − (1 − K(·)) (1/m) and use the empirical distribution function based on K(·). To estimate power of the test with Weibull    distribution as baseline distribution, we generate 10000 samples of sample sizes n = 30, 50 and 100 with parameters (λ, θ) as (1,1), (1,2), (2,1), (2,2), (0.5,1), (0.5,2) and β as 0.5, 0.8, 1.5, 1.8 from the density given in (7). To examine the level attained, we generate 10000 samples of sample sizes n = 30, 50, 100 and 200 with parameters (λ, θ) as (1,1), (1,2), (2,1), (0.5, 1), (0.5, 2) and β = 1. In Tables 5 and 6, respectively, report the level attained and the estimated power of the non-parametric test, with Weibull distribution as baseline distribution. From Table 6 it is observed that the power improves as the sample size increases and as the parameter β moves away from 1. But we need a large sample size to attain 0.05 as level of significance. Similarly, in order to estimate power of the non-parametric test in case of LFR, we generate 10000 samples of sample sizes n = 10, 30 and 50 with parameters (λ 1 , λ 2 ) as (1,1), (1,2), (2,1), (2,2) and β as 0.5, 0.8, 1.5, 1.8 from the density given in (10). To examine at the level attained by the test, we generate 10000 samples of sample sizes n = 30, 50, 200 and 300 with parameters (λ 1 , λ 2 ) as (1,1), (1,2), (2,1), (2,2) and β = 1. We report the level attained and the estimated power of non-parametric test in Table 7 and Table 8 respectively, for LFR baseline distribution. Here also, it is observed that the power improves as the sample size increases and as the parameter β moves away from 1. We note that, to achieve the desired level of significance by the non-parametric test, we need a larger sample size as compare to the LRT. We have used R software for computation purpose.

Illustrations
In present section, we report two illustrations with two data sets namely appendectomy in Australian twins and motors data, to see how the proposed test procedures work.

Appendectomy in Australian twins
The Australian Twins Data (Duffy, Martin, and Mathews (1990)) was derived from a questionnaire mailed to 5967 twin pairs, over the age of 18 years registered with the Australian Twin Registry during the period from November 1980 to March 1982. Subjects were asked if they had undergone appendectomy and at what age the procedure was performed (onset). A total of 3, 808 complete pairs returned the questionnaire. We have taken only the pairs of twins who had undergone appendectomy. There were 339 such pairs. Out of the 339 pairs of twins we discarded 6 pairs of twins for whom the appendectomy was performed at the same age, treating these observation as ties.  We consider a twin pair as a two component parallel system with the age at appendectomy as the time to failure. Generally it is assumed that the twins are associated with each other as far as their health problems are concerned. Our interest is to test whether the age of appendectomy of one individual is associated with the age of the other. We first apply the Wilcoxon-Mann-Whitney type test proposed by Deshpande et al. (2010) to test for the load sharing effect. Here, H 0 is that times of appendectomy of twins are independent and H 1 is that the time of appendectomy of one individual is associated with that of the other. The value of the test statistic is −5.8854 and the p-value = 0.0001, which shows evidence in favor of H 1 . Thus we conclude that the time of appendectomy of one twin has an effect on the time of appendectomy of the other.
Kolmogorov-Smirnov (KS) type test showed that the LFR distribution is appropriate distribution for the component failure time. This test is conservative as we estimated the unknown parameters by the proposed two step procedure. In case of Weibull the two stage estimates of λ, θ and β are 2.0849, 0.0012 and 0.8005 respectively and in case of LFR the two stage estimates of λ 1 , λ 2 and β are 0.0028, 0.0030 and 0.7875 respectively.
In order to investigate whether the load sharing model proposed by Sutar and Naik-Nimbalkar (2016) provides a better fit to the data than the model proposed by Deshpande et al. (2010), we use Akaike Information Criteria (AIC) (Akaike (1974)) and Bayesian Information Criteria (BIC) (Ghosh, Delampady, and Samanta (2007)). We refer to the model proposed by Deshpande et al. where p represents number of parameters of the model, n represents number of data points and L(·) is the maximized value of the likelihood function for the estimated model.
From Table 9, it is clear that LFR model is preferable for the Australian Twins Data. Table  10 shows the AIC and BIC values for the proposed PCRHR and PCHR models with initial component lifetimes as Weibull and as LFR.
From Table 10, we see that, PCRHR model is more appropriate than the PCHR model for the Australian Twins Data. It also shows that the LFR distribution with the PCRHR model gives a better fit.  We also apply the LRT described in Section 3.1 to test the load sharing effect. The value of LRT statistic is 20.6042, which is significant both at 1% and 5% level. Hence we conclude that for a pair, the age of appendectomy of one twin is associated with that of the other twin.

Motors data
These data are taken from Reliability Edge Home (Home (2003)). The data set represents a life test on 18 parallel systems consisting of two electric motors operating continuously. The data consist of the failure times of both the motors along with their identification labels (A and B) for each system. We consider the ordered component failure times for each system. Here, we wish to test whether failures of motors occur independently. Sutar and Naik-Nimbalkar (2016) have shown that the failure of a motor affects the working motor or there exists a load sharing phenomenon. They have also shown that the PCRHR model, proposed by Sutar and Naik-Nimbalkar (2016), is more appropriate than the PCHR model, proposed by Deshpande et al. (2010), for the motors data in case of both Weibull and LFR distributions as the baseline distributions of the components.
We also apply LRT to the same data set. The value of LRT statistic is 6.4038, which is significant both at 1% and 5% level. Hence we may conclude that the failure of a motor affects the working motor. But the value of non-parametric test statistic is 1.7819 which is not significant both at 1% and 5% level, since it is also seen from the simulation study that it reqires larger sample size.

Conclusions
We considered the load sharing model proposed by Sutar and Naik-Nimbalkar (2016). We have proposed a likelihood ratio test and a non-parametric test for tesing the load sharing effect observed in a two component parallel load sharing system. We report the simulation study to check the performance of the proposed test procedures for a two component parallel load sharing system with baseline distribution as Weibull distribution as well as linear failure rate distribution. The simulation study shows that the proposed test procedure is quite satisfactory. To see the practical applicability of the proposed tests, we report the two illustrations.