Two regression methods for estimation of a two-parameter Weibull distribution for reliability of dental materials

Abstract

Objectives

Comparison of estimation of the two-parameter Weibull distribution by two least squares (LS) methods with interchanged axes. Investigation of the influence of plotting positions and sample size. Derivation of 95% confidence intervals (95%CI) for Weibull parameters applicable in the context of LS estimation. Preparation of a free available Excel template for computation of point estimates and 95%CI for Weibull modulus ( m ) and characteristic strength ( s ).

Methods

Monte Carlo simulation covering a wide range of Weibull parameters and sample sizes. Mathematical derivation of formulae for computation of 95%CI according to a Menon-type approach for both m and s . Empirical proof that the practically observed coverage agrees with the nominal one of 95%.

Results

Relative and absolute performance of LS estimators depended on sample size, plotting positions and parameter to be estimated. For most situations they outperformed the corresponding Maximum Likelihood (ML) estimator in terms of bias, while precision was almost the same. Naïve Wald-type 95%CI based on standard errors of LS regression coefficients did not reach targeted coverage. An easy-to-apply alternative based on asymptotic standard errors (Menon 95%CI) resulted in excellent coverage.

Conclusion

Accuracy of the LS methods for Weibull modulus and characteristic strength essentially depend on plotting position and sample size. Large sample sizes ( n ≥ 30) support a credible Weibull parameters estimation. An important complement of the point estimates of Weibull parameters is provided by the Menon 95%CI. A free available Excel template considerably facilitating computation of point and interval estimates of Weibull parameters is provided.

Introduction

Dental restorative ceramics have excellent properties in terms of esthetics and biocompatibility but are susceptible to brittle fracture, a type of failure that is particularly difficult to predict. In general, ceramics are sensitive to defects or flaws within the material, motivating the application of probabilistic concepts and in particular the weakest link (or largest flaw) model. When measuring the flexural strength of a material, the true value may be envisaged as a perfect flaw-free specimen and variations from this true value are due to defects during preparation. Each specimen contains a discrete flaw population and the largest flaw within each test specimen will cause failure if a uniform tension is applied. The distribution of the largest flaws of a sample of test specimens (and therefore the failure strength) follows an extreme value distribution such as Gumbel, Fréchet or Weibull. Among these, the Weibull distribution is most frequently used since it is bounded (lowest possible flexural strength is zero), has a great shape flexibility, can provide accurate failure predictions even with a small number of test specimen, and provides a simple interpretation as a linear regression model (Eqs. (3) and (5) ) .

Weibull distribution

The Weibull distribution allows for prediction of the failure probability at any level of stress providing information about the reliability of a material. The three-parameter Weibull cumulative distribution function (cdf) relates probability of failure G ( x ) = P ( X x ) to stress x :

<SPAN role=presentation tabIndex=0 id=MathJax-Element-1-Frame class=MathJax style="POSITION: relative" data-mathml='G(x)=1−exp−x−s0sm’>G(x)=1exp((xs0s)m)G(x)=1−exp−x−s0sm
G ( x ) = 1 − exp − x − s 0 s m

where s is a scale parameter called characteristic strength, m a shape parameter denoted Weibull modulus and s 0 a threshold representing the minimum stress below a specimen will not break. When the sample size is small ( n < 50, typically the case for dental material investigations) and when there is no physical rationale for a non-zero threshold , it has been suggested to rather consider the two-parameter Weibull distribution with s 0 = 0 . The corresponding cdf and probability density functions (pdf) are:

<SPAN role=presentation tabIndex=0 id=MathJax-Element-2-Frame class=MathJax style="POSITION: relative" data-mathml='G(x)=1−exp−xsm’>G(x)=1exp((xs)m)G(x)=1−exp−xsm
G ( x ) = 1 − exp − x s m
<SPAN role=presentation tabIndex=0 id=MathJax-Element-3-Frame class=MathJax style="POSITION: relative" data-mathml='g(x)=msxsm−1exp−xsm’>g(x)=ms(xs)m1exp((xs)m)g(x)=msxsm−1exp−xsm
g ( x ) = m s x s m − 1 exp − x s m

Unlike for normal distribution, expected values and variance depend on both parameters of the Weibull distribution ( Fig. 1 a and b). Nevertheless, in the context of dental studies the characteristic strength s can be interpreted approximately as the “location” of the distribution, i.e. the stress at which an average specimen breaks. The Weibull modulus m can be seen as an approximate measure of the precision (inverse of the spread) of the distribution. A larger value indicates a narrower distribution and hence a higher precision; a larger fraction of specimen breaks in a smaller range of applied stress. A large Weibull modulus is a highly desirable property for a dental material as it guarantees more uniform performance among different specimens and therefore a higher reliability. The Weibull distribution shows a great shape flexibility ( Fig. 1 a and b) . For 2 < m < 4 the density is fairly symmetric and is well approximated by a normal distribution (for m = 3.4 the best symmetry is attained). A low value of modulus ( m < 1.25) gives a right-skewed, whereas high value of modulus ( m > 10) gives a left-skewed density curve.

Fig. 1
(a and b) Weibull densities for different values of modulus m and characteristic strength s .

Note that there are different parameterization of the Weibull distribution and dependent on the software package other parameters than m and s may be estimated. The R function survreg from package survival , e.g. estimates δ = 1/ m and λ = log( s ), respectively ( Appendix B ).

Monte Carlo simulation

Monte Carlo simulation is a well known approach for performance evaluation of statistical methods . In the context of the two-parameter Weibull distribution it consists of the generation of random numbers (strength data) for a given modulus ( m ) and characteristic strength ( s ).

Example 1 ( Appendix A.3 ): A sample of n = 5 strength values was generated from a Weibull distribution with m and s set to 2 and 10, respectively.

Using Monte Carlo methods, it is certain that the assumption of the underlying Weibull sampling distribution holds and the exact parameter values of this distribution ( m , s ) are known. The advantages of such methodology are manifold. First, strength data from a Weibull distribution with precisely defined m and s can be generated without any necessity for execution of real experiments. Second, parameter values can be freely chosen and, e.g. set to values typical for dental material research. Third, as much data as desirable can be generated and sample size can be varied freely. Fourth, a statistical technique can be applied to a given simulated sample in the same way as to measured data. However, in contrast to measured data, parameter estimates obtained for a simulated sample ( <SPAN role=presentation tabIndex=0 id=MathJax-Element-4-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ,sˆ’>mˆ,sˆmˆ,sˆ
m ˆ , s ˆ
) can be compared to the known true values ( m , s ) to evaluate the performance of the statistical technique in terms of absolute error and bias. The absolute error is defined as <SPAN role=presentation tabIndex=0 id=MathJax-Element-5-Frame class=MathJax style="POSITION: relative" data-mathml='|mˆ−m|’>|mˆm||mˆ−m|
| m ˆ − m |
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-6-Frame class=MathJax style="POSITION: relative" data-mathml='|sˆ−s|’>|sˆs||sˆ−s|
| s ˆ − s |
, respectively. The bias is defined as the difference between the expected value of the estimator (that can be approximated by the mean of <SPAN role=presentation tabIndex=0 id=MathJax-Element-7-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ’>mˆ
m ˆ
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-8-Frame class=MathJax style="POSITION: relative" data-mathml='sˆ’>sˆ
s ˆ
over a large number of samples) and the true values ( m , s ).

<SPAN role=presentation tabIndex=0 id=MathJax-Element-9-Frame class=MathJax style="POSITION: relative" data-mathml='bias(mˆ)=E(mˆ)−m≈mean(mˆi)−mbias(sˆ)=E(sˆ)−s≈mean(sˆi)−s’>bias(mˆ)=E(mˆ)mmean(mˆi)mbias(sˆ)=E(sˆ)smean(sˆi)sbias(mˆ)=E(mˆ)−m≈mean(mˆi)−mbias(sˆ)=E(sˆ)−s≈mean(sˆi)−s
bias ( m ˆ ) = E ( m ˆ ) − m ≈ mean ( m ˆ i ) − m bias ( s ˆ ) = E ( s ˆ ) − s ≈ mean ( s ˆ i ) − s

Error and bias close to 0 indicate that the statistical technique works well. A large absolute error and bias indicate that the statistical technique performs poorly when estimating the Weibull parameters of interest.

Example 2 ( Table A1 ): For YonX (Eq. (3) ) and mean plotting positions (Eq. (7) ), <SPAN role=presentation tabIndex=0 id=MathJax-Element-10-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=1.646′>mˆ=1.646mˆ=1.646
m ˆ = 1.646
with an absolute error of 0.354. In contrast, for XonY (Eq. (5) ) and hazen plotting positions (Eq. (9) ), <SPAN role=presentation tabIndex=0 id=MathJax-Element-11-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=3.948′>mˆ=3.948mˆ=3.948
m ˆ = 3.948
with an absolute error of 1.948.

Estimation of the modulus ( m ) and characteristic strength ( s ) in a two-parameter Weibull distribution

In practice, the least squares (LS) technique for the estimation of the Weibull parameters is preferred due to its simplicity and independence from sophisticated statistical software . As we describe below there are actually two competing methods for LS estimation . We call them here YonX and XonY, respectively. In fact, both approaches can lead to differing estimates of Weibull parameters. What is more, for LS methods the failure probability for each observation ( G ( x i ) = G i , refer to Eq. (1) ) has to be estimated. These failure probability estimates are typically referred to as plotting positions and definition of the plotting positions affects the LS estimates. Below, three definitions for plotting positions will be investigated: mean ranks, median ranks and hazen ranks. Up to now, a clear evaluation of the performance of the two possible LS methods (YonX vs. YonX) and their interplay with plotting positions for varying sample sizes is still missing.

In general, the LS method is based on the linearized form of the Weibull cdf (taking two times the logarithm of Eq. (1) ), allowing for parameter estimation by simple linear regression (ordinary least squares). To fix the notation assume that Y stands for log(−log(1 − G ( x ))) and X for log( x ). Here, log( x ) denotes the logarithm with base e . Let us stress that both LS approaches lead to numerically differing estimates of the Weibull parameters.

Example 3 ( Table A1 , Fig. 2 a ): For YonX (Eq. (3) ), median plotting positions (Eq. (8) ), <SPAN role=presentation tabIndex=0 id=MathJax-Element-12-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=1.911′>mˆ=1.911mˆ=1.911
m ˆ = 1.911
and for XonY (Eq. (5) ), median plotting positions, <SPAN role=presentation tabIndex=0 id=MathJax-Element-13-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=3.431′>mˆ=3.431mˆ=3.431
m ˆ = 3.431
.

Fig. 2
The two different least square linear regression approaches for the data in Appendix A.3 generated for m = 2, s = 10 and n = 5. (a) Comparison of YonX and XonY in the same plot using median ranks for plotting positions. Note that the two regressions do not have identical slopes and intercepts. (b and c) Influence of plotting positions on regression of YonX and XonY.

Least squares (YonX)

This frequently used LS approach resembles the visual approach for Weibull parameters estimation . Its advantage is that the Weibull modulus can be estimated by the slope of the YonX regression directly. The YonX linear regression takes the form Y = b X + a and the correspondence of the intercept ( a ) and slope ( b ) to the functions of m and s can be accomplished by equating the regression parameters with the regression equation obtained from:

<SPAN role=presentation tabIndex=0 id=MathJax-Element-14-Frame class=MathJax style="POSITION: relative" data-mathml='log(−log(1−G(x)))=mlog(x)−mlog(s)→Y=mX−mlog(s)’>log(log(1G(x)))=mlog(x)mlog(s)Y=mXmlog(s)log(−log(1−G(x)))=mlog(x)−mlog(s)→Y=mX−mlog(s)
log ( − log ( 1 − G ( x ) ) ) = m log ( x ) − m log ( s ) → Y = m X − m log ( s )

leading to:

<SPAN role=presentation tabIndex=0 id=MathJax-Element-15-Frame class=MathJax style="POSITION: relative" data-mathml='m=b,s=exp−ab’>m=b,s=exp(ab)m=b,s=exp−ab
m = b , s = exp − a b

Computation of the least squares regression for Y at the ordinate and X at abscissa leads to the estimates <SPAN role=presentation tabIndex=0 id=MathJax-Element-16-Frame class=MathJax style="POSITION: relative" data-mathml='aˆ’>aˆ
a ˆ
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-17-Frame class=MathJax style="POSITION: relative" data-mathml='bˆ’>bˆ
b ˆ
of the regression coefficients and finally to estimates <SPAN role=presentation tabIndex=0 id=MathJax-Element-18-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ’>mˆ
m ˆ
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-19-Frame class=MathJax style="POSITION: relative" data-mathml='sˆ’>sˆ
s ˆ
of the Weibull parameters.

Example 4 ( Table A1 , Fig. 2 a): For YonX (Eq. (3) ), median plotting position (Eq. (8) ), <SPAN role=presentation tabIndex=0 id=MathJax-Element-20-Frame class=MathJax style="POSITION: relative" data-mathml='aˆ=−5.091′>aˆ=5.091aˆ=−5.091
a ˆ = − 5.091
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-21-Frame class=MathJax style="POSITION: relative" data-mathml='bˆ=1.911′>bˆ=1.911bˆ=1.911
b ˆ = 1.911
. Consequently, the Weibull parameter estimates can be obtained: <SPAN role=presentation tabIndex=0 id=MathJax-Element-22-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=bˆ=1.911,sˆ=exp(−aˆ/mˆ)=14.350′>mˆ=bˆ=1.911,sˆ=exp(aˆ/mˆ)=14.350mˆ=bˆ=1.911,sˆ=exp(−aˆ/mˆ)=14.350
m ˆ = b ˆ = 1.911 , s ˆ = exp ( − a ˆ / m ˆ ) = 14.350
. Absolute errors are 0.089 and 4.350, respectively.

Note that in case of YonX the characteristic strength <SPAN role=presentation tabIndex=0 id=MathJax-Element-23-Frame class=MathJax style="POSITION: relative" data-mathml='sˆ’>sˆ
s ˆ
is usually estimated as the uniform stress at which the probability of failure is 0.63 rather than by the formula suggested above <SPAN role=presentation tabIndex=0 id=MathJax-Element-24-Frame class=MathJax style="POSITION: relative" data-mathml='sˆ=exp(−aˆ/mˆ)’>sˆ=exp(aˆ/mˆ)sˆ=exp(−aˆ/mˆ)
s ˆ = exp ( − a ˆ / m ˆ )
.

Least squares (XonY)

The LS regression XonY takes the form X = dY + c with X at ordinate and Y at abscissa. Such a choice is motivated by the fact that X (the actual measurement) varies more than Y and should therefore be the dependent variable in the regression . However, this regression is less straightforward for the practical use since the inverse of the Weibull modulus is estimated by its slope.

<SPAN role=presentation tabIndex=0 id=MathJax-Element-25-Frame class=MathJax style="POSITION: relative" data-mathml='log(x)=1mlog(−log(1−G(x)))+log(s)→X=1mY+log(s)’>log(x)=1mlog(log(1G(x)))+log(s)X=1mY+log(s)log(x)=1mlog(−log(1−G(x)))+log(s)→X=1mY+log(s)
log ( x ) = 1 m log ( − log ( 1 − G ( x ) ) ) + log ( s ) → X = 1 m Y + log ( s )

leading to

<SPAN role=presentation tabIndex=0 id=MathJax-Element-26-Frame class=MathJax style="POSITION: relative" data-mathml='m=1d,s=exp(c)’>m=1d,s=exp(c)m=1d,s=exp(c)
m = 1 d , s = exp ( c )

Computation of the least squares regression for X at the ordinate and Y at abscissa leads to the estimates <SPAN role=presentation tabIndex=0 id=MathJax-Element-27-Frame class=MathJax style="POSITION: relative" data-mathml='cˆ’>cˆ
c ˆ
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-28-Frame class=MathJax style="POSITION: relative" data-mathml='dˆ’>dˆ
d ˆ
of the regression coefficients.

Example 5 ( Table A1 , Fig. 2 a): For XonY (Eq. (5) ), median plotting positions (Eq. (8) ), <SPAN role=presentation tabIndex=0 id=MathJax-Element-29-Frame class=MathJax style="POSITION: relative" data-mathml='cˆ=2.550′>cˆ=2.550cˆ=2.550
c ˆ = 2.550
and <SPAN role=presentation tabIndex=0 id=MathJax-Element-30-Frame class=MathJax style="POSITION: relative" data-mathml='dˆ=0.291′>dˆ=0.291dˆ=0.291
d ˆ = 0.291
. Consequently, the Weibull parameter estimates can be obtained as follows: <SPAN role=presentation tabIndex=0 id=MathJax-Element-31-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ=1/dˆ=3.431′>mˆ=1/dˆ=3.431mˆ=1/dˆ=3.431
m ˆ = 1 / d ˆ = 3.431
, and <SPAN role=presentation tabIndex=0 id=MathJax-Element-32-Frame class=MathJax style="POSITION: relative" data-mathml='sˆ=exp(cˆ)=12.807′>sˆ=exp(cˆ)=12.807sˆ=exp(cˆ)=12.807
s ˆ = exp ( c ˆ ) = 12.807
. Absolute errors of the estimates are 1.431 and 2.807, respectively.

Estimation of G ( x ): plotting positions for LS

A critical point in LS Weibull modeling is the estimation of the failure probability G ( x ). Discrepancies in the LS Weibull statistics found for different software packages can for example be explained by different estimation techniques for G ( x ) . Typically an estimate of the failure probability for each observation (termed <SPAN role=presentation tabIndex=0 id=MathJax-Element-33-Frame class=MathJax style="POSITION: relative" data-mathml='Gˆi’>GˆiGˆi
G ˆ i
) is computed as an order statistic, i.e. a function of the ordered observations ranked from 1 (weakest specimen) to n (strongest specimen). These estimates <SPAN role=presentation tabIndex=0 id=MathJax-Element-34-Frame class=MathJax style="POSITION: relative" data-mathml='Gˆi’>GˆiGˆi
G ˆ i
are sometimes referred to as plotting positions. Among numerous different strategies that have been suggested to estimate failure probabilities, the following three are most frequently used :

<SPAN role=presentation tabIndex=0 id=MathJax-Element-35-Frame class=MathJax style="POSITION: relative" data-mathml='Gˆi=Rin+1,’>Gˆi=Rin+1,Gˆi=Rin+1,
G ˆ i = R i n + 1 ,
<SPAN role=presentation tabIndex=0 id=MathJax-Element-36-Frame class=MathJax style="POSITION: relative" data-mathml='Gˆi=Ri−0.3n+0.4,’>Gˆi=Ri0.3n+0.4,Gˆi=Ri−0.3n+0.4,
G ˆ i = R i − 0.3 n + 0.4 ,
<SPAN role=presentation tabIndex=0 id=MathJax-Element-37-Frame class=MathJax style="POSITION: relative" data-mathml='Gˆi=Ri−0.5n,’>Gˆi=Ri0.5n,Gˆi=Ri−0.5n,
G ˆ i = R i − 0.5 n ,

where R i indicates the rank of the i th observation and n the sample size. These estimates are called mean (Eq. (7) ), median (Eq. (8) ) and hazen ranks (Eq. (9) ). LS linear regressions of YonX and XonY using the different plotting positions are shown in Fig. 2 b and c. Mean ranks have been suggested early and are widely used. Median ranks have been claimed to be most accurate , whereas hazen ranks have been shown to lead to the least biased estimates for medium to large sample sizes ( n > 20) . In particular, mean ranks have been used in combination with regression of YonX and median ranks in combination with XonY . Bias caused by plotting positions in the context of the XonY estimation is also discussed in .

Example 6 ( Table A1 , Fig. 2 b): Errors and bias can also be caused by the plotting positions definition within one LS estimation technique. For example for YonX regression: <SPAN role=presentation tabIndex=0 id=MathJax-Element-38-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ’>mˆ
m ˆ
(mean ranks) = 1.646, <SPAN role=presentation tabIndex=0 id=MathJax-Element-39-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ’>mˆ
m ˆ
(median ranks) = 1.911, <SPAN role=presentation tabIndex=0 id=MathJax-Element-40-Frame class=MathJax style="POSITION: relative" data-mathml='mˆ’>mˆ
m ˆ
(hazen ranks) = 2.164 can differ considerably. Absolute errors of the estimates are 0.354, 0.089 and 0.164, respectively.

Maximum likelihood (ML)

ML is an alternative estimation technique. It has attractive mathematical properties such as consistency, asymptotic normality and efficiency and allows a simple construction of confidence intervals . ML estimates are calculated by maximizing the likelihood of m and s given the data.

<SPAN role=presentation tabIndex=0 id=MathJax-Element-41-Frame class=MathJax style="POSITION: relative" data-mathml='argmaxm,sL(m,s|x1,…,xn)=argmaxm,s∏i=1ng(xi)=argmaxm,s∏i=1nmsxism−1exp−xism’>argmaxm,sL(m,sx1,,xn)=argmaxm,sni=1g(xi)=argmaxm,sni=1ms(xis)m1exp((xis)m)argmaxm,sL(m,s|x1,…,xn)=argmaxm,s∏i=1ng(xi)=argmaxm,s∏i=1nmsxism−1exp−xism
arg max m , s L ( m , s | x 1 , … , x n ) = arg max m , s ∏ i = 1 n g ( x i ) = arg max m , s ∏ i = 1 n m s x i s m − 1 exp − x i s m
Only gold members can continue reading. Log In or Register to continue

Nov 23, 2017 | Posted by in Dental Materials | Comments Off on Two regression methods for estimation of a two-parameter Weibull distribution for reliability of dental materials
Premium Wordpress Themes by UFO Themes