Why didn’t you ask before running all those tests?

The title is a Statistician’s Lament. “Design of Experiments (DoE) is the design of any task that aims to describe or explain the variation of information under conditions that are hypothesized to reflect the variation.” [Wikipedia] Are you using DoE to design reliability tests? What do PH, GMDH, and |D|-optimality have to do with design of DoE of reliability tests?

The W79 nuclear artillery shell needed additional safing to prevent criticality when stored or when loaded into M115 or M110 howitzer. [https://en.wikipedia.org/wiki/W79_Artillery-Fired_Atomic_Projectile/] The safing mechanism was to be removable material inside the shell’s pit; arming removed the material allowing implosion criticality while preserving shell safety [Washington Post May 23, 1990]. A chemist did 350 tests of alternative materials, quantities, and removal methods, without conclusion about shell safety. PH (proportional hazards), GMDH (Global Method of Data Handling, www.gmdh.net), plus sixty more |D|-optimal DoE tests assured mechanical safing reliability [NIST, Boon].

Steve Fisher of Marvin Windows asked for meta-reliability analysis of sealant test failure times data [Sugiyama]. Failure times varied from minutes to 9 years, many irrelevant to window lives. The sealant tests were not well-designed, so I combined the PH model and the GMDH method, and used least-squares to fit the PH-GMDH model at failure-time order statistics. I used the model to estimate sealant reliability for a 10-year warranty.

PH-GMDH Model of Failure Rate Functions

The proportional hazards model in biostatistics, clinical trials, and credible reliability prediction assumes the failure rate functions are λ(t)=λ_o(t)*Exp[zβ], where z is a vector of “concomitant” variables representing test factors, λ_o(t) is the “common” failure rate function (when z=0), β is a vector of regression coefficients, and zβ = SUM[z(i)*β(i), i = 1,2,…,k] [Bobrowski et al.]. PH models are convenient, because they combines a nonparametric failure rate function with parametric proportions Exp[zβ]. The nonparametric failure rate function doesn’t require unwarranted assumptions, and the regression coefficient magnitudes β(i) indicate relative importance of the test factors (if the z-vector components are standardized). Fan and Li show how to select explanatory variables in PH and “frailty” models. 

Why not use polynomial instead of linear regression zβ: for better fits, exploratory data analysis, quantifying effects and interactions of sealant and load factors, and DoE? The GMDH method for polynomial regression is a form of artificial intelligence [Farlow], especially if it is not apparent how z-variables influence reliability. The GMDH method doesn’t require unwarranted assumptions, physics-of-failure, and it might be explainable to management. Consider using PH-GMDH models for life tests, optimal Design of Experiments, and reliability predictions and analyses.

PH and GMDH applications are not new. “Credible Reliability Prediction” [George] uses a PH model, λ(t; new)  λ(t; old)*MTBF(old)/MTBF(new), because, during design of a new product, all that is known is MTBF(new), the field failure rate function(s), and MTBF(old(s)) of comparable older parts and products. Generations of product failure rate functions are proportional, because parts, processes, environments, and customers change little. “Frailty” PH-models in survival analysis use alternatives to linear zβ regression. 

GMDH uses the Ivakhnenko polynomial [Farlow], 

IP(z) = a+Σb_ix_i+ΣΣc_ijx_ix_j+ΣΣΣd_ijkx_ix_jx_k+…,

with IP(z) as a function of independent variables x, where the z-vector denotes products of x variables. The sums are over all sample and all indexes. (IP(z) denotes that polynomial in the rest of this article. It’s customary in DoE and regression to describe y as a function of x, and in PH to denote the variable combinations as z-vector. PH-GMDH uses IP(z) instead of linear zβ in a Proportional Hazards reliability model. 

Sealant Test Data

Sealant test factors included: M50 the secant tensile modulus at 50% extension, Tmax the maximum tensile strength, and Emax (%) the elongation at Tmax and hours-to-failures in steady load. There were 24 samples of three different cure systems and several combinations of filler systems labeled A, B, C, D, E, F, G. The test data are in the Test Engineering and Management article by George and in the SugiPHaz.xlsx workbook.

Sugiyama graphed histograms, computed sample means, and linearly extrapolated mean lives to low loads. Except for extrapolation, this is OK for the tested sealants and loads, within ranges tested. NIST, ASTM C24 and C510, SAE, Dow [Wolf, Sugiyama] specify strength test methods and data collection, but not reliability test methods. 

Reliability is the probability of survival to specified ages under specified conditions, P[Life > t]. Reliability affects warranty, spares, and uncertainty costs, and contributes to revenue. Profit should motivate interest in reliability, but it doesn’t seem to justify the sealant test design.

• Some test were at extremely high loads? The company’s 10-year warranty  87,660 hours. 
• Sample coefficient of variation (CV=standard deviation/mean) varied from 18.9% to 198%, for shear loads. High CV is undesirable for product lives. Sample skew and kurtosis were both positive and negative.
• Two test loads do not justify linear extrapolation; at least three loads are needed to test linearity of mean lives. Some mean lives aren’t linear.
• Some test samples were at same loads in both shear and tensile? Strengths differ in shear and tensile. Reliabilities in shear and tensile differ even more.
• Reliability people care about P[Life > warranty] and P[Life > useful life], not just mean life.

PH-GMDH Method

The objective is to model sealants’ reliability for different loads and strengths, using the test data and an explanatory subset of measurable variables. The model could also be used to design more reliability tests for other conditions and sealants; i.e., apply the model to other sealants’ stresses (cyclic, temperature, etc.). 

Treat observed lives as sample time-to-failure order statistics: k = 1, 2,…, 12. Estimate each of the 12 sample’s reliability function percentiles R(t_k) as {11/24, 13/24,…,123/24} at order statistics, and estimate the corresponding actuarial failure rate function estimates as λ(t_k) (R(t_k1)R(t_k))/(t_k1). Not every sample had 12 distinct failure times, so count repeated failure times as successive order statistics. Some tests ended before failure, so use the Nelson-Aalen estimator for right-censored samples [Wikipedia].

The GMDH method determines polynomial powers and estimates the coefficients. GMDH uses a step-wise regression on combinations of independent variables to explain variations in the data. It tries different polynomial models until goodness-of-fit is not improved. It tries pairwise polynomials b₁z_i+b₂z_j+b₃z_i²+b₄z_j²+b₅z_iz_j for all pairs (z_i, z_j) and the remaining variables. Then it replaces original variables with the best-fit (z_i, z_j) polynomial. At each iteration, estimate the GMDH polynomial coefficients and common failure rate function by minimizing the goodness-of-fit objective:
SUM(observed-expected)²/expected = SUM[λ(t_k; z)λ_o(t_k)*Exp[IP(z)])²/(λ_o(t_k)*Exp[IP(z)])].

The sum is over k = 1,2,…,12 and over all samples. Iterate until fit deteriorates, run out of variables, or AIC (Akaike Information Criterion) is too big [Farlow]. Figure 1 shows λ_o(t_k) for alternative z-variables. Winner λ_o(t_k)Exp{IP(z)] uses z = {Load, M50, Load₂, M50², Load*M50).

The SugiPHaz.xlsx workbook includes:

• “Table2Data” includes some exploratory data analyses: empirical reliability functions, correlations, CVs, Weibull fits, regresses Weibull parameters on factors, RPH plots
• “TensilePH” does regression and PH-GMDH analyses
• “ShearPH” sets up data for multi-factor PH-GMDH
• “ShearPHEmax”, “ShearPHM50”, and “ShearPHTmax” do PH-GMDH analyses
• “TenShear” models both tensile and shear data
• “Summary” collects SSEs, graphs l_o(t_k) and R_o(t) estimates, fits some distributions to estimates, and generates nonparametric order statistics for any l_o(t_k)*Exp[IP(z)]

The Mathematica CoxModelFit[] function uses different methods to fit rhw function form of the PH proportionality factor such as Exp[IP(z)]. Mathematica package MmaGMDH is a “program for the implementation of the traditional GMDH algorithm,” which could provide the function form.  Combining CoxModelFit[] and MmaGMDH could improve on my ad hoc workbook method.

Using Results

The PH-GMDH model of reliability is P[Life > t] = EXP[λ_o(s)*Exp[IP(z)ds], where the variable of integration is from 0 to t. Collect PH-GMDH model parameter estimates, graph λ_o(t_k) and R_o(t), optionally fit parametric reliability functions, and generate order statistics for any failure rate λ_o(t_k)*Exp[IP(z)]. 

Commonly assumed parametric reliability functions don’t fit sealant bimodal distributions very well (figure 2). (Stochastic ordering of failure rate functions yields parametric reliability functions, sometimes. [Kochar]) From the SugiPHaz.xlsx workbook, retrieve IP(z) coefficients from TensilePH tables 4 and 5 or ShearPHM50 tables 3 and 4. Use “Summary” spreadsheet table 3 or 5, λ(t_k; z) = λ_ο(t_k)Exp(IP(z)) for your z-vector. Generate nonparametric order-statistic t_k-values by solving for t_k in 

λ_ο(t_k)*Exp(IP(z)) = (R(t_k1)R(t_k))/((t_kt_k1)*R(t_k1)). 

Design PH-GMDH Reliability Tests

|D|-optimal design maximizes |X’X|, the determinant of “information” matrix X’X, where X is the “design” matrix of independent variable values; (e.g., variables from tables 1 and 2 or their GMDH polynomial variables z. |D|-optimal design selects subsets of variable values to maximize information from a specified number of tests and the GMDH model. [NIST Engineering Statistics Handbook section 5.5.2.1, Boon] For example, to add one more tensile test to either samples B or C, |D|-optimally, either test one more B-sample at 0.625 or at 0.417 MPa. The former was derived using Load and Tmax, and the latter was derived using Load, Tmax, Load², Tmax², and Load*Tmax. The differences between |X’X|-values were negligible, limited by computational accuracy. 

|D|-optimal design minimizes a “generalized” variance of parameter estimates, whereas we’d like to minimize expected errors E[R(t)Exp[-Exp[IP(z)]*λ_o(u)du]], (integrate u from 0 to t), between reliability R(t) and the PH-GMDH model for some t- and z-values. Use |D|-optimality until you get more specific about expected error costs. E.g.; w = warranty, find z-values so that Exp{IP(z)]ln(1/R(w))/λ_o(u)du (integrate from 0 to w).

Use PH-GMDH model to estimate loads and reliability

Steve Fisher of Marvin WIndows would like sealant G to survive a 10-year warranty; what shear and tensile loads could be tolerated? The sealant G shear test data indicate that it is unlikely to survive a year, using the simple PH failure rate model λ_o(t)*Exp[β*Load]. Friend said that sealant F was similar but stronger. Sealant F has 50% probability of surviving 10 years, if shear load is 0.11 MPa.  

The sealants B and C tensile test data were used to estimate the tensile load for a 10-year warranty. The Ivakhnenko polynomial used Load and M50 with linear Tmax and Emax. Ranges of samples B and C lives were so large that there was no chance of high 10-year reliability, but 0.025 MPa tensile load would allow ~50% of sample F to survive 10 years. 

Technically, window survival requires surviving both tensile and shear loads, as well as other insults. The Sugiyama test data indicates 25% probability of surviving both 0.11 MPa shear and 0.025 MPa tensile loads for 10 years, assuming independence. Marvin Windows’ field reliability is much better. 

Recommendations

Use the Ivakhnenko polynomial IP(z) in the proportional hazards model of failure rate function λ_o(t)Exp[IP(z)], where d IP(z)/dz is the rate of change of the (logarithm of the) proportionality factor. Equate failure rate functions at order statistics, for modeling reliability of different materials with different loads. Combine cyclic loads, temperature, etc. Use physically relevant material strength parameters; e.g. Emax seems irrelevant for tensile load reliability. Combined shear and tensile model doesn’t fit as well as separate: tensor physics is not reliability. Test within the operating ranges of loads and lives. If not, use some reasonable, physically justified model, or else try the GMDH method, and design tests optimally to verify reliability. Send test information or data and ask questions, or request workbook from pstlarry@yahoo.com. 

References

Bobrowski, S., M. Döring, U. Jensen, and W. Schinköthe, “Reliability Prediction Using the Cox Proportional Hazards Model,” 56-th International Scientific Symposium, Ilmenau University of Technology, Sept. 2011

Boon, Jacob E., “Generating Exact D-Optimal Designs for Polynomial Models,” SpringSim, Vol. 2, 2007

Fan, Jianqing and Runze Li, “Variable Selection for Cox’s Proportional Hazards Model and Frailty Model,” The Annals of Statistics, Vol. 30, No. 1, pp. 74–99, 2002  

Farlow, Stanley J., “A FORTRAN Program for the GMDH Algorithm,” Chapter 15 of Self-Organizing Methods in Modeling,CRC Press, 1984  

Farlow, Stanley J., “The GMDH Algorithm of Ivakhnenko,” The American Statistician, Vol. 35, No. 4, pp. 210-215, Nov. 1981

George, L. L., “Credible Reliability Prediction,” 2^nd Edition, See CREDRP2020.pdf in  https://sites.google.com/site/fieldreliability/credible-reliability-prediction/.

George, L. L., “What Were You Thinking? PH-GMDH Sealant Reliability Model,” Test Engineering and Management,Vol. 78, No. 5, Oct.-Nov., 2016

Kochar, Subhash C., “Stochastic Comparisons of Order Statistics and Spacings: A Review,” ISRN Probability and Statistics, vol. 2012, Article ID 839473, 2012 

NIST, Engineering Statistics Handbook, http://www.itl.nist.gov/div898/handbook/pri/section5/pri521.htm

Sugiyama, S., “Creep Resistance of Structural Glazing Silicone Sealants,” Dow Corning Toray Silicone Co., Chiba, Japan

Wolf, A. T, “Durability Testing of Sealants,” Sealant Technology Conference, Oxford Brooks Univ., Oct. 2004