Fitting a distribution to my timeseries: two R-packages, two contrary results

Question

typical R-User here, applying a bunch of packages to my data, hoping for a convincing result although I understand only half of what I do at most (and then getting rejected by the reviewers, surprise). But bear with me, I am here to learn.

I have a time-series of monthly precipitation data, for which I'd like to find a suitable distribution. The Cullen and Frey graph from the fitdistrplus package indicates that a lognormal distribution should outperform a gamma distribution:

and the results of fitting a gamma, Weibull and lognormal distribution seem to say the same:

Goodness-of-fit statistics
                               Weibull  lognormal      gamma
Kolmogorov-Smirnov statistic 0.0263855 0.03245692 0.07115254
Cramer-von Mises statistic   0.1427163 0.11374430 1.13149886
Anderson-Darling statistic   1.0934641 0.62867455 6.73258465

Goodness-of-fit criteria
                                Weibull lognormal    gamma
Akaike's Information Criterion 7845.461  7842.426 7954.480
Bayesian Information Criterion 7854.940  7851.905 7963.959

This is of particular interest to me, as in my field of research a gamma distribution is often used for precipitation data, without considering that there might be better fitting distributions.

I was quite satisfied with my results until I stumbled upon the GAMLSS package, which offers an all-in-one solution for finding a suiting distribution to your data. So I applied the easy to use `fitDist' function to see what's the outcome:

> fit <- fitDist(df$prcp, k = round(log(length(df$prcp))), type = "realplus", trace = F , try.gamlss = TRUE)
There were 11 warnings (use warnings() to see them)
> summary(fit) 
*******************************************************************
Family:  c("BCCG", "Box-Cox-Cole-Green") 

Call:  gamlssML(formula = y, family = DIST[i], data = sys.parent()) 

Fitting method: "nlminb" 


Coefficient(s):
            Estimate  Std. Error   t value   Pr(>|t|)    
eta.mu    45.1741563   0.9560784  47.24943 < 2.22e-16 ***
eta.sigma -0.5522153   0.0265535 -20.79636 < 2.22e-16 ***
eta.nu     0.4250894   0.0439835   9.66474 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   842 
Global Deviance:     7827.23 
            AIC:     7833.23 
            SBC:     7847.45 

First of all, the function offers me a distribution I've never heard of, but that's fair since (you might've already guessed it) I am quite new to this statistical topic it makes sense that I don't know "all" distributions. But, as I investigated the results of the other distributions, it got quite interesting (values are the BIC):

> fit$fits
    BCCG    BCCGo       GA       GG      BCT     BCTo     BCPE    BCPEo     WEI2      WEI     WEI3      GB2      GIG   exGAUS      LNO 
7848.230 7848.230 7852.426 7853.005 7854.763 7854.763 7854.989 7854.989 7855.461 7855.461 7855.461 7857.338 7859.426 7871.134 7964.480 
  LOGNO2    LOGNO       IG      EXP PARETO2o  PARETO2       GP   IGAMMA 
7964.480 7964.480 8260.617 8290.105 8297.106 8297.106 8297.112 8521.495 

Right behind the mysterious BCCG-distributions follows the GA (gamma) distribution, quite closely followed by Weibull, and then, lognormal ranked even lower.

So here are my questions:

1. What's the reason for these contrary results? Is it due to different fitting methods?
2. What can I further do to find the most fitting distribution?

I know in general that all models are wrong and that it probably won't make a huge difference whether I pick the result of fitdistrplus or GAMLSS for my further analysis (calculation of SPI). But still, I am wondering about the comparability of research:

3. In a publication, I would probably state the package that I used for the fitting and present some goodness-of-fit values (such as the BIC e.g.). But, obviously, there's more to that and the results depend on the package used. But in the publications that I investigated, there's rarely more information offered. Doesn't that hinder the comparability of research?

Answer 1

With the outlier included in the data, the maximum likelihood estimation of the parameters for a particular distribution may be converging to a local maximum, rather than the global maximum.

You could try different starting values for the parameters, e.g. fit each distribution separately and include argument sigma.start=0.5 in the gamlss() function.

Without the outlier, the fits should be more stable, and so less chance of a local maximum. Does the contradiction you found disappear when the outlier is removed?

Answer 2

Another idea is to use weights to weight out the outlier. Then use the resulting fitted model (say m1) as starting values for the distribution parameters, when fitting the model (say m2) with all the data points (i.e. without weights).

This is done by using argument, e.g. start.from=m1, in the function

m2 <- gamlss( ,start.from=m1)