Inconsistency between R and SAS for MLE on Weibull

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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.