Solved – fitting a log-logistic distribution

distributionsestimationlogisticr

I am conducting survival analysis research and I am trying to decide which distribution (Weibull, lognormal and loglogistic) best fits my data.

One can "Fit of univariate distributions to non-censored data" using the fitdist() from the fitdistrplus package in R.

referring to the examples from the help page, we can fit hypothetical distributions by the following:

fitW <- fitdist(serving, "weibull")
fitg <- fitdist(serving, "gamma")
fitln <- fitdist(serving, "lnorm")

etc.

But one cannot fit an assumed log-logistic distribution with the package.

From research, I have found this paper (see pg 5-6) which attempts to implement a secondary package, actuar which has the capability to fit a loglogis distribution.

Here we find the code:

fendo.ll <- fitdist(ATV, "llogis", start = list(shape = 1, scale = 500))

which then states at the bottom of the page

It can thus help to find correct initial values for the distribution parameters in non trivial cases, by iterative calls if necessary (see the reference manual for examples (Delignette-Muller et al., 2014))

So can one explain where these values come from?

Best Answer

I am not sure I understand your problem, but this can be done directly with fitdistr. See the following code:

library(fitdistrplus)
library(actuar)

test <- actuar::rllogis(500, 3.3, 1)
fitdist(test, "llogis")
Fitting of the distribution ' llogis ' by maximum likelihood 
Parameters:
      estimate Std. Error
shape 3.374404 0.12585423
scale 0.998631 0.02296057

The fitdistr function do not have a set of predefined distribution models, it just assumes that pdqr... functions corresponding to basename llogis exist somewhere on the search path.

Related Solutions

Solved – confirm or validate underlining distribution association with survival analysis

The use of log-likelihoods are problematic because not all the survival models are nested within one another. But you can use something like AIC or BIC, which can be obtained using extractAIC(). For example:

> extractAIC(survival.weibull)
[1]    5.000 2284.504
> extractAIC(survival.logn)
[1]    5.000 2315.142
> extractAIC(survival.logl) 
[1]    5.000 2297.171

This would suggest, based on AIC, that the Weibull model fits best. In terms of visualizing the best fit, you could plot all three parametric fits atop a non-parameteric Kaplan-Meier curve to visualize the fit of the parametric forms to the underlying data. The means to do that is in your first question.

Solved – What (discrete) distribution might fit best here

Let me open by saying that searching for and identifying a model that appears to describe the data is not where I would normally begin my modelling; it begins with the information whuber and I were seeking on how the model will be used, with what is being measured - and further, to consider any associated theoretical or expert practical understanding of its properties, the performance of models on similar data sets, and so on. Considerations like these are essential to the choice of an appropriate model.

Nevertheless, as requested, I will attempt to describe some steps one could take in examining data like this and seeing if it could suggest a model.

If you're going to use the resulting distributional model for inference (hypothesis tests, or confidence intervals, or prediction intervals for example), there's a problem that you chose the model based on your data so the nominal properties of the inference will no longer hold. (If you really can't use the characteristics of the variable and subject area knowledge, having no other basis than the sample itself, then it's better to pull off a random subset and use that for model identification.)

Given the heavy tail, one option is to look at the way the log of the empirical survivor function relates to $x$. My first plots were to plot $\log(\hat{S}(x))$ against $\log(x)$ (not shown; linear at the top but then somewhat curved near the middle) and then against $\sqrt x$. The appearance of these then suggested to me that a fair first model might be a zeta distribution. However, they can also suggest a number of other possibilities (such as perhaps directly modelling the seeming linear relationship in the first plot below; it might provide a somewhat better fit, at the expense of being less familiar)

So I then looked at the empirical pmf, plotting $\log(\hat{p}(x))$ against $\log{x}$ as a check on the potential for fitting a zeta function.

That indeed suggests a fairly straight line relationship, and it looks like a parameter in the ballpark of 2 or perhaps a bit more would be a reasonable first model -- note that the most precisely estimated proportions are the first few, and fitting those accurately will matter more for a lot of applications -- but not all. However, it's not a perfect fit - with a large sample size like you have, I believe you would easily reject a pure zeta model (or indeed any other simple model).

(NB the parameter value here has not been optimized; it's just a reasonable guess at a suitable value)

So if I had to suggest a simple model on the basis of only looking at the data, the zeta distribution would be an obvious candidate; if an initial guess at the parameter value is needed, I'd start at 2.

On fitting power laws, such as the zeta distribution, also see Cosma Shalizi's essay So you think you have a power law? which debunks some of the enthusiasm for power laws (a lot of things will look more or less power-law like) and the paper it discusses, which has a lot of useful information.

Best Answer

Related Solutions

Solved – confirm or validate underlining distribution association with survival analysis

Solved – What (discrete) distribution might fit best here

Related Question