Inconsistency between R and SAS for MLE on Weibull

I have the following data Y which I want to get an MLE estimate for the parameters using a Weibull distribution in R.

1468, 1872, 475, 1372, 3830, 1849, 978, 1389, 909, 701, 1227, 962, 1781, 580, 584, 2675, 841, 1544, 452, 955, 556, 1737, 747, 1565, 1331, 1188, 2649, 1800, 2718, 808, 1138, 909, 1359, 846, 1334, 1397, 719, 1715, 681, 2002, 994, 2543, 1564, 1717, 1106, 1859

If I try and run fitdistr(Y, "weibull") I get a warning:

fit = fitdistr(Y, "weibull")
Warning message:
In densfun(x, parm[1], parm[2], ...) : NaNs produced
> warnings(fit)
Warning message:
In densfun(x, parm[1], parm[2], ...) : NaNs produced Error in
at(list(...), file, sep, fill, labels, append) :
argument 2 (type 'list') cannot be handled by 'cat'

But it still gives me a MLE. However, the value is different from the result SAS gives.

output from R:

shape          scale
2.1103684   1537.2344072
(0.2245888)  (112.1596367)

output from SAS (using proc lifereg):

Weibull Scale 1550.559
Weibull Shape 2.1195

What's causing the discrepancy and are there any preferred packages/functions for computing simple estimates of MLEs for distributions over MASS and fitdistr?

• It is impossible to get a "precise" point estimate, that's the thinking behind using confidence intervals. In Bayesian case you don't bother with point estimates and treat both the data and model parameters as random variables with certain distributions. – Tim Nov 16 '14 at 9:02
• The difference between the two sets of parameters is much smaller than the bracketed estimates of the errors in these estimates. So they are close. – Henry Nov 16 '14 at 16:06
• I have two suggestions: i) As the parameters of the weibull distributions have to be greater than or equal to zero, you should add the option lower = 0 to the call of fitdistr (so the error messages of the optim function disappear). ii) Supply better starting values for shape and scale, say 10 and 150. Using this setup with fitdistr, the parameter estimates were: shape = $2.1262279$, scale = $1563.1042469$. These are very close to @Yves estimates. – COOLSerdash Nov 16 '14 at 16:34
• Stata program weibullfit (SSC) gives 2.126228 and 1563.095, i.e. very close to results of @COOLSerdash. – Nick Cox Nov 16 '14 at 17:13
• I have found fitdistr to be very fragile and you will get much better results writing a quick negative logliklihood function in R and passing it to optim or nloptr. – Avraham Jan 29 '15 at 5:12

First remind that the fitdistr function (from the MASS package) is a very general function that can work with nearly any distribution. The warnings come from non-allowed parameter values (e.g. negative scale or shape) met during the optimisation unconstrained by default.

It seems a good idea here to try a specific MLE for the Weibull distribution. A quite well-known fact is that the ML estimation of the two-parameter Weibull can be rely on a concentration of the log-likelihood, leading to an easier one-dimensional optimisation. Moreover, the concentrated log-likelihood is concave, so there is a unique ML estimate.

The problem here is that the log-likelihood is quite flat near the optimum, so different optimisations lead to different results as reported by @Glen_b. Moreover, the data scaling is prone to numerical problems. After rescaling, similar results are obtained with or without concentration. A general practical finding about MLE is that using poorly scaled data can be enough to ruin the estimation.

> library(Renext)            ## for concentrated log-lik
> try(fweibull(Y))           ## error (numerical pb with information matrix)
> fit <- fweibull(Y / 1000)  ## works
> ## set parameters and logLik back to original scale
> fit$est * c(1, 1000) shape scale 2.126225 1563.094460 > fit$sd * c(1, 1000)
shape       scale
0.2444308 114.1293266

> fit$loglik - length(Y) * log(1000) [1] -362.2237 > library(MASS) > ## set parameters and logLik back to original scale > fit2 <- fitdistr(Y / 1000, "weibull") > fit2$est * c(1, 1000)
shape       scale
2.126231 1563.095165

> fit2$sd * c(1, 1000) shape scale 0.2288605 114.9071653 > fit2$loglik - length(Y) * log(1000)
[1] -362.2237

An optimization function shouldn't be expected to give identical answers to a similar function in a different package - or even to the same function with different options.

I tried a variety of different optimizers and starting places in fitdistr. They generally gave really similar results, of which the SAS and the result you got with fitdistr in R were typical.

I have included one of those fits, using a different optimizer in fitdistr and non-default starting point. In terms of the resulting fit, all three are essentially indistinguishable (and your two results are more alike than the third):

I don't think anything is amiss.

The warning should not be ignored, but investigated as far as possible, but sometimes errors (or in this case, warnings) can be generated without indicating that there's any convergence problem. You should try to figure out what caused that. Trying different starting points and optimizers (and plotting the resulting fits) should indicate if there's much of an issue.

[Ideally, you should plot the function in a 3D plot (or its contours in a 2D plot) near the identified optimum, which will help identify a number of potential problems.]

With the Weibull, one thing you can do is use the survreg function in the survival package, which will fit a Weibull as its default model. Its two parameters are related to the usual Weibull ones (this is described in the help on survreg). You just want a constant-mean model:

> survreg(Surv(Y)~1)
Call:
survreg(formula = Surv(Y) ~ 1)

Coefficients:
(Intercept)
7.354423

Scale= 0.4703164

Loglik(model)= -362.2   Loglik(intercept only)= -362.2
n= 46
> exp(7.354423)   #  exponentiate the Intercept
[1] 1563.095
> 1/0.4703164     #  take inverse of the Scale
[1] 2.126228

summary(survreg) will give standard errors on the scale it uses, but if you take say a 95% CI and transform the endpoints, they can be used as a CI for the transformed parameters.

While the SAS output is better than the R output, the unpleasant fact is that both perform rather poorly. To see this, note that the gradient at the reported solution should disappear to 0 ... whereas for both the R and SAS results, this is not the case.

In particular, let $$X \sim Weibull(b,c)$$ with pdf $$f(x)$$:

(source: tri.org.au)

I am going to activate mathStatica's SuperLog function:

(source: tri.org.au)

Then, the exact symbolic log-likelihood for $$\theta = (b,c)$$ is given by:

(source: tri.org.au)

Replacing $$(x_1, \dots, x_n)$$ by the $$n=46$$ data values yields the exact observed log- likelihood:

(source: tri.org.au)

For the R solution, and the SAS solution, the gradient vector, calculated at each reported solution, is:

(source: tri.org.au)

At the optimal solution, the gradient should disappear to 0. The SAS solution is better than the R solution, but both are poor. The solution reported by Yves does much better:

(source: tri.org.au)

... but can still be easily improved upon.

The Hessian matrix (at the solution) ... and the eigenvalues of the Hessian ... should also be calculated to ensure that the observed log-likelihood is concave in the neighbourhood.

• wolfies, where does the gamma distribution enter the picture? Am I missing something? – COOLSerdash Nov 16 '14 at 16:41
• @COOLSerdash Oops - just a typo. The density entered is the Weibull and all is fine ... I have fixed the name typo. Thanks for picking that up. – wolfies Nov 16 '14 at 16:52