Paper de referencia: Delignette-Muller and Dutang 2015 Fitdistrplus: An R package for fitting distributions 64: 1-34.

Describe como ajustar datos a distribuciones. Para ello utiliza el paquete fitdistrplus que permite ajustar distribuciones univariadas de diferentes tipos de datos, y además permite diferentes métodos de estimación (maximum likelihood, moment matching, maximum goodness-of-fit estimation, etc.). Además se pueden comparar diferentes ajustes entre si. Esta herramienta es muy útil en ecotoxicología y en evaluación del riesgo ambiental.

Manual de fitdistrplus.

Instalamos el paquete fitdistrplus y sus dependencias install.packages("fitdistrplus")

FITTING DISTRIBUTIONS TO CONTINUOUS NON-CENSORED DATA

Para ello se utiliza el ejemplo groundbeef. Datos de consumo de hamburguesas de carne (gr) en niños menores de 5 años en Francia.

library("fitdistrplus")
## Loading required package: MASS
## Loading required package: survival
library("MASS")
library("survival")
data("groundbeef", package="fitdistrplus")
str(groundbeef)
## 'data.frame':    254 obs. of  1 variable:
##  $ serving: num  30 10 20 24 20 24 40 20 50 30 ...

Lo primero es observar la distribución empírica (si no se sabe cuál debería seguir teóricamente). Con las siguientes dos figuras se tiene el histograma de frecuencias y la distribución acumulada. Para ello:

plotdist(groundbeef$serving, histo=TRUE, demp=TRUE)

Ahora algunos estadísticos nos pueden ayudar a decidirnos por el tipo de distribución, entre ellos la kurtosis y skewness. Si el valor de skewness es diferente de cero indica que carece de simetría la distribución empírica, y la kurtosis cuantifica el peso de los extremos, ya que un valor de 3 es la kurtosis de una distribución normal. Para ello se emplea descdist:

descdist(groundbeef$serving)

## summary statistics
## ------
## min:  10   max:  200 
## median:  79 
## mean:  73.64567 
## estimated sd:  35.88487 
## estimated skewness:  0.7352745 
## estimated kurtosis:  3.551384

FIT OF DISTRIBUTIONS BY MAXIMUM LIKELIHOOD ESTIMATION

Para ajustar nuestros datos a alguna de las distribuciones teóricas de los datos se utiliza la función fitdist. Cualquier distribución tiene que tener definidos tres parámetros (d, p, q), que son density, distributin y quantile functions. Para saber más de las distribuciones ?Distributions.

El argumento dist indica para la función fitdist la distribución a la que se quiere ajustar los datos. Por ejemplo norm para una distribución normal.

La función fitdistnos da la siguiente información:

-Los parámetros estimados.

-Errores estándard.

-Loglikelihood.

-Akaike and Bayesian criteria (AIC y BIC).

-Matriz de correlación entre los parámetros estimados

Un ejemplo con la distribución weibull. Más información sobre eta distribución en la web de minitab.

fw<-fitdist(groundbeef$serving, "weibull")
summary(fw)
## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##        estimate Std. Error
## shape  2.185885  0.1045755
## scale 83.347679  2.5268626
## Loglikelihood:  -1255.225   AIC:  2514.449   BIC:  2521.524 
## Correlation matrix:
##          shape    scale
## shape 1.000000 0.321821
## scale 0.321821 1.000000
plot(fw)

Se pueden realizar las figuras por separado de la siguiente forma:

denscomp(fw)

cdfcomp(fw)

qqcomp(fw)

ppcomp(fw)

Ahora toca comparar varias distribuciones teóricas para los datos de este ejemplo. En este caso se hará para weibull, lognormaly gamma, de la siguiente manera:

fg<-fitdist(groundbeef$serving, "gamma")
fln<-fitdist(groundbeef$serving, "lnorm")
par(mfrow=c(2,2))
plot.legend<-c("Weibull","lognormal","gamma")
denscomp(list(fw,fln,fg), legendtext=plot.legend)
qqcomp(list(fw,fln,fg), legendtext=plot.legend)
cdfcomp(list(fw,fln,fg), legendtext=plot.legend)
ppcomp(list(fw,fln,fg), legendtext=plot.legend)

Según la interpretación de los autores del manuscrito “The density plot and the CDF plot may be considered as the basic classical goodness-of-fit plots. The two other plots are complementary and can be very informative in some cases. The Q-Q plot emphasizes the lack-of-fit at the distribution tails while the P-P plot emphasizes the lack-of-fit at the distribution center. In the present example, none of the three fitted distributions correctly describes the center of the distribution, but the Weibull and gamma distributions could be preferred for their better description of the right tail of the empirical distribution, especially if this tail is important in the use of the fitted distribution, as it is in the context of food risk assessment”.

EJEMPLO CON DATOS ECOTOXICOLÓGICOS

Para ello se utilizará el ejemplo de endosulfan, que es un insecticida organoclorado. Todos los valores se expresan en microgramos por litro. Para más infrmación sobre el endosulfan.

Los datos del ejemplo provienen de Hose y Van den Brink 2004) en un trabajo de laboratorio con especies australianas. Se trata de una SSD (Species Sensitivity Distribution). En estas distribuciones de sensibilidades se suelen calcular las HC5 (Hazardous Concentrations 5%), que es la concentración de tóxico que protege al 95% de las especies del ecosistema. Lo primero que se debe hacer es ver como es la distribución de estos datos:

data("endosulfan", package="fitdistrplus")
str(endosulfan)
## 'data.frame':    104 obs. of  3 variables:
##  $ ATV       : num  0.6 2.8 182.2 0.8 478 ...
##  $ Australian: Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 1 ...
##  $ group     : Factor w/ 3 levels "Arthropods","Fish",..: 1 1 1 1 1 1 1 1 1 1 ...
plotdist(endosulfan$ATV, histo=TRUE, demp=TRUE)

descdist(endosulfan$ATV)

## summary statistics
## ------
## min:  0.1   max:  39891.8 
## median:  4.9 
## mean:  1522.39 
## estimated sd:  5720.916 
## estimated skewness:  5.076456 
## estimated kurtosis:  30.728

Para estos datos la distribución lognormalo loglogistic no se ajustan muy bien (son las más habituales para SSDs). Para conocr esto hacemos lo siguiente (hay que instalar el paquete actuar):

data("endosulfan", package = "fitdistrplus")
ATV <- endosulfan$ATV
fendo.ln <- fitdist(ATV, "lnorm")
library("actuar")
## 
## Attaching package: 'actuar'
## The following object is masked from 'package:grDevices':
## 
##     cm
fendo.ll <- fitdist(ATV, "llogis", start = list(shape = 1, scale = 500))
fendo.P <- fitdist(ATV, "pareto", start = list(shape = 1, scale = 500))
fendo.B <- fitdist(ATV, "burr", start = list(shape1 = 0.3, shape2 = 1, rate = 1))
cdfcomp(list(fendo.ln, fendo.ll, fendo.P, fendo.B), xlogscale = TRUE,
ylogscale = TRUE, legendtext = c("lognormal", "loglogistic", "Pareto",
"Burr"))

La parte derecha de la distribución queda bien descrita por todas las distribuciones elegidas, pero la izquieda no. Ésta es la importante, ya que hay que calcular la HC5, de las cuatro seleccionadas la distribución Burr es la que mejor se ajusta. Para calcular la HC5 se hace del siguente modo:

quantile(fendo.B, probs = 0.05)
## Estimated quantiles for each specified probability (non-censored data)
##             p=0.05
## estimate 0.2939259

Y el cuantil 0.05 teórico:

quantile(ATV, probs = 0.05)
##  5% 
## 0.2

Ahora toca determinar estadísticamente la bondad del ajuste de los datos a los diferentes modelos. Para ello hay tres estadísticos Cramer-von Mises, Kolmogorov-Smirnovy Anderson-Darling (D´Agostino y stephens 1986). Para utilizarlos empleamos la función gofstat:

gofstat(list(fendo.ln, fendo.ll, fendo.P, fendo.B),fitnames = c("lnorm", "llogis", "Pareto", "Burr"))
## Goodness-of-fit statistics
##                                  lnorm    llogis     Pareto       Burr
## Kolmogorov-Smirnov statistic 0.1672498 0.1195888 0.08488002 0.06154925
## Cramer-von Mises statistic   0.6373593 0.3827449 0.13926498 0.06803071
## Anderson-Darling statistic   3.4721179 2.8315975 0.89206283 0.52393018
## 
## Goodness-of-fit criteria
##                                   lnorm   llogis   Pareto     Burr
## Akaike's Information Criterion 1068.810 1069.246 1048.112 1045.731
## Bayesian Information Criterion 1074.099 1074.535 1053.400 1053.664

El estadístico de Anderson-Darlingda más peso a los extremos, y para ERA es el preferido. En todo caso hay que quedarse con los valores más bajos, tanto de estadístico como de Aikake y Bayesian. En este ejemplo parece que los ajustes de Burr y Pareto son los mejores.

Ahora hay que calcular el grado de incertidumbre de los modelos. Se puede hacer por bootstraps paramétricos o no paramétricos con la función boodist. Lo haremos para la distribución Burr del ejemplo de endosulfan:

bendo.B <- bootdist(fendo.B, niter = 1001)
summary(bendo.B)
## Parametric bootstrap medians and 95% percentile CI 
##           Median      2.5%     97.5%
## shape1 0.2018955 0.1005319 0.3772913
## shape2 1.5832546 1.0374096 2.8267319
## rate   1.4702089 0.6582622 2.7495060
plot(bendo.B)

Si queremos intervalos para quantiles concretos la forma de hacerlo es la siguiente:

quantile(bendo.B, probs = 0.05)
## (original) estimated quantiles for each specified probability (non-censored data)
##             p=0.05
## estimate 0.2939259
## Median of bootstrap estimates
##             p=0.05
## estimate 0.3025125
## 
## two-sided 95 % CI of each quantile
##           p=0.05
## 2.5 %  0.1708425
## 97.5 % 0.5070410

ADVANCED TOPICS

En el artículo se analizan varias aproximaciones para ajustar modelos, este apartado es interesante para distribuciones complicadas, incluso se pueden customizar el algoritmo óptimo para la distribución de los datos.

También se pueden aplicar a datos semi-cuantitativos.

REFERENCIAS

-D’Agostino RB, Stephens MA (1986). Goodness-of-Fit Techniques. 1st edition. Dekker.

-Hose GC, Van den Brink PJ (2004). Confirming the Species-Sensitivity Distribution Concept for Endosulfan Using Laboratory, Mesocosm, and Field Data. Archives of Environmental Contamination and Toxicology, 47(4), 511-520.