1
$\begingroup$

I am trying to reproduce a plot from a paper. The plot is a Weibull probability plot and is used to illustrate why a Weibull distribution is appropriate for modelling failure times.

The plot is shown below.Weibull Probability Plot

In order to reproduce this image. I assumed that the failure times followed a Weibull distribution. I obtained the ML estimates and corresponding 95% CIs for the parameters (the values obtained matched the values obtained by the authors).

I then calculated the probability of failing for times between 0 and 70 using the Weibull CDF with the ML estimates, and with the lower and upper points of the CIs.

I then plotted the actual data points (proportion of units failing over time) and the probability of failing using the three pairs of parameters described above. My plot is shown below.

My plot

My plot does not look like the plot the authors produced. Have I misunderstood what is meant my a probability plot?

I understand that the axes chosen by the author are different to my axis, but if you focus on the x-axis between 0 and 70, all of their data points seem to be inside the CI. This is not the case for my plot.

Edit: Below are the ML estimates and corresponding 95% CIs for the parameters of the failure time distribution.

\begin{equation} \begin{aligned} Parameter& \hspace{2cm}MLE& \hspace{2cm} Lower&\hspace{2cm} Upper& \\ \mu& \hspace{2cm} 6.741120& \hspace{2cm} 6.0660978& \hspace{2cm} 7.416142&\\ \sigma& \hspace{2cm} 0.9029375& \hspace{2cm} 0.7228623& \hspace{2cm} 1.127872& \end{aligned} \end{equation}

Here $T \sim Weibull(\mu, \sigma)$, where the cdf is given by

\begin{equation} F(t, \mu, \sigma) = 1 - \exp\bigg(1-\exp\bigg(\frac{log(t) - \mu}{\sigma}\bigg)\bigg) \end{equation}

The red line represents a Weibull cdf evaluated at values between 0 and 70 with $\mu = 6.0660978$ and $\sigma = 0.7228623$. The green line a Weibull cdf evaluated at values between 0 and 70 with $\mu = 7.416142$ and $\sigma = 1.127872$.

The paper I am referring to is

Further Edit: To clarify: The code and data used to create the (second) plot is shown below. I (currently) do not know how to make the axes as shown in the first plot, but I will attempt to do this.

#############################################################################################
#Phi_{sev}
psev=function(z)
{
  1-exp(-exp(z))
}
#phi_{sev}
dsev=function(z)
{
  exp(z-exp(z))
}
#############################################################################################
#Creating functions which will evaluate the Weibull pdf and cdf 
FF = function(x, mu, sigma){
  psev((log(x) - mu)/sigma)
}

ff = function(x, mu, sigma){
  (1/(x*sigma))*dsev((log(x) - mu)/sigma)
}

#############################################################################################
#MLE estimates
MLE = c(6.7411199, 0.9029375)
MLE_lower = c(6.0660978, 0.7228623)
MLE_upper = c(7.416142, 1.127872)

tt = seq(from = 0.1, to = 70, length.out = 10000)

Weibull_MLE = FF(tt, MLE[1], MLE[2])
Weibull_lower = FF(tt, MLE_lower[1], MLE_lower[2])
Weibull_upper = FF(tt, MLE_upper[1], MLE_upper[2])

#############################################################################################
#Data
failure.time = c(0.3522425,  1.8514547,  2.1562701,  3.4413818,  4.7055079,  4.8174972,  5.0486417,  5.4441475,  5.7243062,
                 5.9805326,  6.4711694,  7.0957963,  7.7831483,  8.7863524, 10.7897027, 10.9652542, 11.7657082, 12.5167305,
                 15.1997829, 16.0759416, 16.3112283, 16.5817208, 17.1930128, 17.9957668, 18.6290973, 18.6963592, 20.2800036,
                 21.1053799, 22.1641800, 22.1807384, 22.2426650, 22.6613307, 23.4989884, 24.8859476, 24.9962985, 25.5614063,
                 25.7941020, 26.3491067, 27.0350730, 27.5901930, 27.7006820, 28.1287065, 28.2125510, 28.8704583, 32.1762028,
                 32.5667568, 32.8155859, 33.7885900, 36.3435663, 37.0979577, 38.6961335, 39.6225545, 39.7377064, 40.7381252,
                 41.6626021, 42.6420496, 45.0285338, 45.1072121, 45.4133994, 46.4729561, 47.0328201, 48.2121433, 49.1853123,
                 50.1805967, 51.1500872, 53.2659645, 53.4411921, 56.2414894, 64.5524150)

prop.fail = 1:69/1800

#############################################################################################
#Creating plot

plot(failure.time, prop.fail, xlab = "Weeks after Installation", ylab = "Fraction Failing", ylim = c(0, 0.05))
lines(tt, Weibull_MLE)
lines(tt, Weibull_lower, col = 2)
lines(tt, Weibull_upper, col = 3)
legend("topleft", legend=c("MLE", "Lower", "Upper"),
       col=c("black", "red", "green"), lty = c(1,1,1))
$\endgroup$
  • $\begingroup$ Compare the labels on the vertical axes. $\endgroup$ – whuber Sep 11 '19 at 17:58
  • $\begingroup$ I fail to see that all of the data points are inside the CI. The green and red lines cross around the 40th week and only data points up to the 30th week seem to be within them. $\endgroup$ – polettix Sep 11 '19 at 17:59
  • $\begingroup$ Edited. I mean all of their data points are inside the CI, but mine are not. That is why I am assuming they have done something different to what I understood to be a probability plot. $\endgroup$ – JLee Sep 11 '19 at 18:03
  • 2
    $\begingroup$ Note that the lower limit of your 95% CI is above the upper limit of your CI past ~38. I'm guessing there are some issues with your calculations at a minimum. Without seeing all of that, it's hard to know why this went wrong. $\endgroup$ – gung - Reinstate Monica Sep 11 '19 at 18:28
  • 2
    $\begingroup$ 1. Ignore the confidence interval on the plot for the moment and let us begin by getting the basic plot right. (i) Note that the probability plot at the top has log time on the horizontal axis, whereas your plot appears to be linear in time. (ii) as whuber noted already, the vertical axis on the top plot is also not linear; yours looks like it is linear. 2. Your approach to the 95% interval may contain some misunderstandings - it would help if you were to provide a more detailed explanation of what you did and how that gives a (presumably pointwise) interval for the QQ plot-points. $\endgroup$ – Glen_b -Reinstate Monica Sep 12 '19 at 2:21