La estadística descriptiva tiene más que ver con describir una muestra de manera cualitativa (gráfica) o cuantitativa (numérica) que con inferir propiedades acerca de la población de la cual proviene esa muestra (estadística inferencial).
La estadística descriptiva está fuertemente ligada al análisis exploratorio de datos (EDA for Exploratory Data Analysis) y al análisis inicial de datos (IDA por Initial Data Analysis). El primero focaliza en explorar los datos en busca de nuevas hipótesis, las cuales pueden terminar en nuevos muestreos y experimentos, mientras el segundo se focaliza en asegurar la calidad de los datos, chequear las asunciones y realizar las transformaciones necesarias para testear la hipótesis que teníamos en mente a la hora de recolectar los datos.
El análisis inicial de datos no es necesario solamente para el testeo de hipótesis de la estadística inferencial, sino también como un paso previo para el aprendizaje automático (ML por Machine Learning) y forma parte importante de las primeras fases de la minería de datos (Data mining).
Se denomina muestra a un subconjunto de datos tomados o seleccionados de una población estadística mediante un proceso de muestreo determinado. Cada una de las unidades muestrales suele llamarse observación, y es posible medir variables aleatorias para cada una de ellas.
Las muestras pueden ser:
En el muestreo aleatorio simple se selecciona un número k de unidades muestrales de manera aleatoria, teniendo cada elemento de la población la misma probabilidad de ser seleccionado. El muestreo puede ser con o sin reposición (with or without replacement). Si el muestreo se efectúa con reposición, cada elemento de la población puede ser seleccionado más de una vez. En el muestreo aleatorio sin reposición, donde los elementos no son devueltos a la población y no pueden ser elegidos más de una vez, la probabilidad de sacar un determinado elemento cambia con la extracción del anterior (no son independientes). Sin embargo el muestreo aleatorio simple sin reposición satisface intercambiabilidad, es que cualquier orden de los elementos extraídos es igualmente probable. Si el tamaño de la población es mucho más grande que el tamaño de la muestra, el muestreo aleatorio simple sin reposición se aproxima a un muestreo simple con reposición, dada la baja probabilidad de elegir un mismo elemento dos veces.
La manera más sencilla de hacer un muestreo aleatorio simple con o sin reposición en Julia, es usando la función sample de StatsBase. Ésta es similar a la función sample de R. Ambas toman una lista de valores, el tamaño de la muestra a tomar (size) y se les debe indicar si queremos que sea con o sin reemplazo o reposición (replace). Una diferencia entre ambas funciones es que mientras R por defecto hace el muestreo sin reemplazo, Julia lo hace con reemplazo si no le indicamos lo contrario.
In [1]:
individuos = 1:10
Out[1]:
In [2]:
using StatsBase
sample(individuos, 5, replace=true)
Out[2]:
In [3]:
using RCall
R"sample($individuos, size=5, replace=TRUE)"
Out[3]:
El muestreo sistemático consiste en ordenar los elementos de una población según alguna variable de interés, para luego tomar n unidades muestrales equiespaciadas. El primer elemento debe ser seleccionado al azar, quedando los otros determinados en relación a este. Una ventaja de este tipo de muestreo es que permite muestrear una variable de intereses en todo su rango. Pero debe tenerse cuidado si la variable muestra alguna característica periódica, dado que puede generarse una muestra sesgada de la población. Además debe tenerse en cuenta que no se verá la variación entre dos elementos contiguos en la lista, dado que nunca son seleccionados a la vez.
Por ejemplo, si queremos hacer un muestreo sistemático de las plantas del conjunto de datos iris según largo de sus pétalos:
In [5]:
using RDatasets
iris = dataset("datasets","iris")
head(iris)
Out[5]:
In [6]:
sort!(iris, cols=:PetalLength) # Ordena la variable de interés
head(iris)
Out[6]:
In [7]:
N = nrow(iris) # Número de elementos en mi población
Out[7]:
In [8]:
n = 7 # Número de muestras
Out[8]:
In [9]:
@assert n != 0 "El tamaño de la muestra no puede ser 0"
@assert n <= N "La muestra no puede ser superior en tamaño a la población"
k = div(N,n) # Intervalo de muestreo (sampling interval or skip)
Out[9]:
In [10]:
primero = sample(1:k) # Elijo al azar (con igual probabilidad) un elemento de 1 a k para comenzar
Out[10]:
In [15]:
# Selecciono los elementos restante a k pasos del primero
índices = Int[ primero + ((muestra-1) * k) for muestra in 1:n ]
Out[15]:
In [16]:
iris[índices,:]
Out[16]:
En el muestreo estratificado, se estratifica la población antes de tomar las muestras. En este proceso se divide a los miembros de la población en estratos, grupos o subpoblaciones homogéneas. Estos estratos deben ser mutuamente excluyentes, dado que un miembro de la población no podrá pertenecer a más de una subpoblación. Todos los miembros de la población deben pertenecer a un estrato determinado, no pueden quedar miembros sin clasificar (estratificación exhaustiva). Una vez que se estratifico la población, se realizar un muestreo aleatorio simple o sistemático dentro de cada estrato.
Existen tres posibles estrategias:
Por ejemplo, si consideramos que todas las filas del conjunto de datos iris es nuestra población:
In [17]:
μₛ = mean(iris[:PetalLength])
Out[17]:
In [20]:
estratos = by(iris, :Species) do df
DataFrame(
μₕ = mean(df[:PetalLength]),
Nₕ = nrow(df)
)
end
Out[20]:
In [21]:
N = nrow(iris)
Out[21]:
In [22]:
(1/N) * sum(estratos[:μₕ] .* estratos[:Nₕ])
Out[22]:
Los estadísticos de resumen (Summary statistics) describen de manera cuantitativa la distribución de una muestra. Normalmente se obtiene un conjunto de estadísticos que describen cada variable aleatoria/dimensión de los datos de manera independiente (si excluimos medidas de dependencia como las correlaciones).
Existen estadísticos de:
Estos estadísticos pueden ser robustos o no, llamándose robustos a los estadísticos menos afectados por valores atípicos (outliers). Estos estadísticos describen la muestra de una manera más robusta cuando su distribución se aleja de una distribución normal (por ejemplo, si la distribución es asimétrica).
Existen diversos estadísticos de tendencia central, siendo la media el más popular de ellos a pesar de no ser un estadístico robusto. La medidas de ubicación robusta más popular es la mediana. La moda es el único estadístico de tendencia central para datos nominales, pero es difícil de estimar correctamente para variables continuos.
In [23]:
for (nombre, valores) in eachcol(iris)
println(nombre)
if eltype(valores) <: Real
println(" Media : ", mean(valores))
println(" Mediana : ", median(valores))
println(" Moda : ", mode(valores))
else
println(" Moda : ", mode(valores))
end
end
De los estadísticos de dispersión, la desviación estándar es el más popular acompañando a la media. No debe confundirse con el error estadístico o estándar, que en realidad habla de la dispersión de las medias muestrales, y no de la variable de interés. La desviación estándar no es un estimador robusto. Alternativas más robustas son el rango entre cuartiles y la Median Absolute Deviation (MAD).
In [24]:
for (nombre, valores) in eachcol(iris)
if eltype(valores) <: Real
println(nombre)
println(" Varianza : ", var(valores))
println(" Desviación estándar : ", std(valores))
println(" Rango intercuartílico : ", iqr(valores))
println(" Median Absolute Deviation : ", mad(valores))
end
end
La kurtosis de la distribución normal es 3, sin embargo suele usarse la $kurtosis - 3$ o excess kurtosis para comparar contra la distribución normal. Distribuciones para las cuales es menos probable observar outliers que en la distribución normal, dan valores de kurtosis negativa. Mientras que valores positivos de kurtosis se observan para distribuciones de probabilidad para los cuales se obtienen outliers o valores extremos con más frecuencia (es decir, poseen colas más pesadas).
In [47]:
using Distributions
using Plots
plt = plot()
for distribución in [ Uniform(-2,2), Normal(), Laplace() ]
long_name = string(distribución)
short_name = match(r"\.(\w+)\{",long_name)[1] # Obtiene Normal de Distributions.Normal{Float...
println(long_name)
println(" Kurtosis: ", kurtosis(distribución))
plot!(plt, -3:0.01:3, x -> pdf(distribución, x), label=short_name)
end
plt
Out[47]:
In [25]:
for (nombre, valores) in eachcol(iris)
if eltype(valores) <: Real
println(nombre)
println(" Skewness : ", skewness(valores))
println(" Kurtosis : ", kurtosis(valores))
end
end
In [53]:
using StatPlots
pyplot(size=(600,300))
Out[53]:
In [54]:
columnas = [:SepalLength, :SepalWidth, :PetalLength, :PetalWidth]
Out[54]:
Los Histogramas permiten tener una visión de cómo sería la forma de una distribución de densidad para una variable aleatoria continua. Se construyen dividiendo la variable en grupos (bins) y contando el número de observaciones dentro de cada uno (representada por la altura de la barra).
El siguiente paso de complejidad que podríamos dar para tener una mejor estimación de la función de densidad de probabilidad, es utilizar Averaged Shifted Histograms (ASH) o el estimador por núcleo (KDE por Kernel Density Estimator). En Julia, si la biblioteca KernelDensity está instalada, Plots la utilizará en su función density. R también posee una función density que puede observarse usado plot(density(....
KDE tiene dos parámetros importantes, uno es la función kernel a utilizar que deber ser una distribución de probabilidad, por defecto se utiliza la distribución Normal. El otro parámetro es el ancho de banda a utilizar.
In [55]:
plot(
histogram(iris, columnas, alpha=0.5, bins=0:0.2:8, legend=true),
density(iris, columnas, alpha=0.5, legend=true, fill=0)
)
Out[55]:
Los Diagramas de Dispersión (scatter plots) utilizan las coordenadas cartesianas para mostrar cómo se distribuyen dos variables en un espacio bi dimensional. Es posible representar más dimensiones utilizado diferentes formas, tamaños y/o colores.
In [61]:
plot(
scatter(iris, :SepalLength, :SepalWidth, group=:Species),
scatter(iris, :PetalLength, :PetalWidth, group=:Species, legend=false)
)
Out[61]:
Otra manera de observar la distribución conjunta de dos variables continuas es haciendo uso de los histogramas bivariados o bidimensionales (2D o bivariate histograms). Los grupos (bins) se establecen para las dos variables, definiendo rectángulos en un espacio bidimensional. Normalmente se utiliza un código de color para indicar la cantidad de valores en cada grupo.
In [64]:
histogram2d(iris, :PetalLength, :PetalWidth, size=(400,300), bins=10)
Out[64]:
Esta representación se asemeja a la de un Mapa de calor (heatmap), pero estos últimos generalmente poseen variable discretas o categóricas en sus ejes (en lugar de la discretización de una variable continua, aunque ésto también es posible). Los heat maps son útiles para mostrar los valores de una matriz, y representan a una variable continua utilizando un gradiente de colores.
In [66]:
medias = by(iris, :Species) do df
DataFrame(
SepalLength = mean(df[:SepalLength]),
SepalWidth = mean(df[:SepalWidth]),
PetalLength = mean(df[:PetalLength]),
PetalWidth = mean(df[:PetalWidth])
)
end
Out[66]:
In [81]:
heatmap(string.(columnas), string.(medias[:Species]), convert(Matrix{Float64}, medias[:,columnas]),
ratio=:equal, yflip=true)
Out[81]:
Muchas veces es útil pasar un DataFrame de un formato wide a uno long. En esta conversión se genera una nueva variable categórica con el nombre de las columnas como niveles y se acumula sus valores en una nueva variable. En Julia es posible usar melt y stack para cambiar la forma del DataFrame de wide to long o unstack para ir en el otro sentido. Esta operación es particularmente útil en R cuando se trabaja con ggplot2, y requiere de la biblioteca reshape2.
In [82]:
long = melt(iris, :Species)
head(long)
Out[82]:
Los Diagramas de Cajas (boxplots) son otra forma de visualizar datos acerca de una distribución. Son interesantes porque permiten observar estadísticos robustos de tendencia central (mediana) y de dispersión (rango entre cuartiles).
In [94]:
boxplot(string.(long[:variable]), long[:value], notch=true, legend=false)
Out[94]:
The notches (if requested) extend to +/- 1.58 IQR/sqrt(n). This seems to be based on the same calculations as the formula with 1.57 in Chambers et al (1983, p. 62), given in McGill et al (1978, p. 16). They are based on asymptotic normality of the median and roughly equal sample sizes for the two medians being compared, and are said to be rather insensitive to the underlying distributions of the samples. The idea appears to be to give roughly a 95% confidence interval for the difference in two medians. R boxplot's help
In [85]:
# install.packages("ggplot2") # Para installar ggplot2 en R
In [95]:
R"""
library(ggplot2)
ggplot($long, aes(x=as.character(variable), y=value)) + geom_boxplot()
"""
Out[95]:
Gadfly es una biblioteca similar a ggplot2 pero para Julia.
In [97]:
using Gadfly
In [98]:
Gadfly.plot(long, x=:variable, y=:value, Geom.boxplot)
Out[98]:
Mientras los density plots permiten observar una PDF inferida a partir de los datos, la función ecdf de la biblioteca StatsBase de Julia,
permite observar la función de densidad de probabilidad acumulada empírica de una muestra (Empirical Cumulative Distribution Function ECDF), una aproximación a la CDF.
In [99]:
using StatsBase
In [100]:
ECDF = ecdf(iris[:SepalLength]) # retorna una función
Out[100]:
In [101]:
plot(ECDF, extrema(iris[:SepalLength])..., label="ECDF")
Out[101]:
Es posible, para los datos de una muestra, ajustar una distribución de densidad de probabilidad (distribution fitting). Muchos de los estadísticos determinados durante la exploración de datos pueden son los estimadores de máxima verosimilitud o maximum-likelihood estimator (MLE) de algunas distribuciones. Por ejemplo, la media es el estimador de máxima verosimilitud de la media poblacional $\mu$ de una distribución normal. Tanto los estadísticos cuantitativos, como el análisis gráfico realizado en la etapa de exploración (histogramas, density plots, boxplots y ecdf) pueden ayudar a decidir cuál familia de funciones de densidad de probabilidad debemos ajustar. La biblioteca Distributions de Julia posee un método fit que permite estimar los parámetros de una distribución a partir de los datos, en general usando máxima verosimilitud.
In [102]:
using Distributions
variable = iris[:SepalWidth]
mean(variable)
Out[102]:
In [103]:
std(variable) # std usa n-1
Out[103]:
In [104]:
ajuste = fit(Normal, variable) # std usa n
Out[104]:
In [106]:
PDF = histogram(variable, normed=true)
plot!(PDF, x -> pdf(ajuste,x), extrema(variable)...)
CDF = plot(ecdf(variable), extrema(variable)...)
plot!(CDF, x -> cdf(ajuste,x), extrema(variable)...)
plot(PDF, CDF, legend=false)
Out[106]:
Como parte del análisis inicial de datos, muchas veces debemos comprobar si los valores de una variable en nuestra muestra se distribuyen de manera normal. Existen diversas manera de testear normalidad, las más populares son:
Los gráficos QQ (qq plots por quantile) permiten comparar los cuantiles de dos distribuciones. Cuando en uno de los ejes se colocan los cuantiles de una distribución normal canónica, hablamos de un gráfico de probabilidad normalidad. Si la distribución de los datos es aproximadamente normal, el gráfico de normalidad mostrará los puntos equivalente a los datos distribuidos sobre una línea.
In [111]:
Gadfly.plot(x=ajuste, y=variable, Stat.qq, Geom.point)
Out[111]:
Existen diversos test de hipótesis para testear normalidad (Ghasemi & Zahediasl 2012), estos tienen como hipótesis nula que los datos provienen de una cierta distribución (en particular de una distribución normal). Por lo tanto, p values altos son indicativos de que no hay evidencia suficiente para rechazar que los datos provienen de una distribución normal. Las dos pruebas de normalidad recomendadas son la de Anderson-Darling y la prueba de Shapiro-Wilk, ambas poseen más poder que Kolmogorov-Smirnov para esa tarea.
Ghasemi, Asghar, and Saleh Zahediasl. "Normality tests for statistical analysis: a guide for non-statisticians." International journal of endocrinology and metabolism 10.2 (2012): 486-489.
In [113]:
using HypothesisTests
In [114]:
OneSampleADTest(variable, Normal())
Out[114]:
In [115]:
R"shapiro.test($variable)"
Out[115]:
Se denomina teorización post hoc a la generación de hipótesis sugeridas por el conjunto de datos observado, sin testear esta hipótesis en nuevos datos. Hacerlo puede resultar en aceptar hipótesis incorrectas, que sólo son válidas en el presente conjunto de datos, dado que es más probable aceptar una hipótesis testeada en el conjunto de datos en el cual se genero.
Es necesario testear estas nuevas hipotesis en una nueva muestra de la población. Sin embargo, en muchas casos eso puede ser imposible, por ejemplo al analizar un fenómeno natural finito. Este problema de disponer de un conjunto limitado de datos, haciendo difícil o imposible la recolección de nuevos datos para la fase de confirmación fue denominado por John Tukey como uncomfortable science.
Un caso no tan obvio de uncomfortable science o de teorización post hoc puede surgir bioinformática cuando se toman todos los elementos disponibles de una base de datos para hacer los análisis exploratorios sin dejar datos suficientes para testear las hipótesis que surgen de ese análisis.
Otro proceso que puede generar este tipo de problemas es el conocido como data fishing, el cual consiste en ir testeando hipótesis sobre un conjunto de datos hasta encontrar una un valor estadísticamente significativo en algún test de hipótesis. El término data fishing posee una connotación negativa, dado que representa el comportamiento poco ético de buscar, seleccionar e informar sólo los resultados positivos sin el debido control, muchas veces como si las hipótesis hubiera sido anterior al análisis (Leung 2011).
Buscar patrones en los datos, muchas veces aplicando test estadísticos, no es incorrecto. Lo incorrecto es utilizar el mismo conjunto de datos (in-sample data) para soportar esa hipótesis. Las alternativas posibles para no caer en este problema son:
Leung, Kwok. "Presenting post hoc hypotheses as a priori: Ethical and theoretical issues." Management and Organization Review 7.3 (2011): 471-479.