La estadística inferencial trata de inferir información acerca de la distribución de probabilidad de una población, partiendo de la información de una muestra. Se diferencia de la estadística descriptiva en que en lugar de simplemente describir la muestra, trata de describir algo acerca de la población. Uno puede realizar diferentes estimaciones acerca de la población, la cual normalmente se conduce haciendo uso de la estadística frecuentista:
Una prueba de hipótesis es análoga a una demostración matemática por contradicción (en particular es un tipo de reductio ad absurdum), aunque no debe considerarse en ningún caso una demostración matemática. Ambos métodos, parten de considerar una declaración verdadera, en el caso del test de hipótesis, la hipótesis nula ó H₀. En ambos métodos se trata de llegar a una contradicción, dado que una implicación solo puede ser falsa si partiendo de una premisa verdadera llegamos a una conclusión falsa (tabla de verdad):
| p | q | p → q |
|---|---|---|
| true | true | true |
| true | false | false |
| false | true | true |
| false | false | true |
Acá es donde ambos métodos divergen, dado que la demostración matemática lo hace siguiendo reglas lógicas determinadas y procedimientos matemáticos permitidos, llegando a un resultado determinado no azaroso. Finalmente, si se llega a una contradicción, se concluye que la declaración original es falsa, aceptando la negación de dicha declaración. En el caso de la prueba de hipótesis, rechazar H₀ implica aceptar la hipótesis alternativa Hₐ. Sin embargo, no rechazar H₀ no implica aceptarla, decimos simplemente que no hay evidencia suficiente para rechazar la hipótesis nula. En cualquier caso, una prueba de hipótesis no demuestra una hipótesis, sólo la soporta con un cierto grado de confianza. Acumular evidencia a favor de una hipótesis (razonamiento inductivo) no es suficiente para demostrarla, la evidencia acumulada solo puede soportar dicha hipótesis.
Las pruebas de hipótesis fueron propuestas para solucionar el problema de la inducción e incluir el análisis estadístico dentro de un modelo hipotético deductivo (Kaplan 2019, capítulo 13).
Existen dos tipos de errores que podemos cometer (Neyman & Pearson, 1933):
| H₀ es cierta | Hₐ es cierta | |
|---|---|---|
| Se escogió Hₐ | Error de tipo I | Decisión correcta |
| Se escogió H₀ | Decisión correcta | Error de tipo II |
El riesgo o probabilidad de caer en un error del tipo I se conoce como nivel de significación estadistica, normalmente indicado con la letra griega $\alpha$. La probabilidad de caer en un error del tipo II normalmente se anota con la letra griega $\beta$. Siendo $1-\beta$ la potencia o poder del test.
Se habla de errores de tipo III cuando se da la respuesta correcta (rechazar una hipótesis nula) para una pregunta incorrecta o errónea.
El valor P (p value) es la probabilidad de observar un determinado valor (o uno más extremo) dado que H₀ es correcta. En ningún caso debe interpretarse como una probabilidad de la hipótesis nula. Un valor P bajo indica que los datos son inconsistentes con la hipótesis nula, permitiéndonos rechazar la hipótesis nula ante esa observación. Para decidir que es un valor bajo, uno puede recurrir al valor de significación α.
La prueba de T de Student para una muestra permite probar si la media de una muestra es diferente de una media $\mu_0$. Algo que tenemos que tener en cuenta, es que no podemos testear que la media poblacional sea un determinado valor. Ésto se debe a que no podemos saber cuál es la media poblacional sin analizar toda la población. Es decir, no podemos saber dónde se encuentra localizada la distribución T Student correspondiente a la media poblacional. Para determinar probabilidades necesitamos un parámetro de ubicación exacto. Es por esto que la hipótesis nula será la ubicación de la media poblacional en un lugar determinado. De hecho, las hipótesis nulas se caracterizan normalmente por tener igualdades (uno supone una distribución con parámetros determinados).
Si para una muestra de tamaño $n$ se cumplen las supuestos del test:
Entonces el estadístico $t$ sigue una distribución T de Student con $n - 1$ grados de libertad.
Supongamos una variable aleatoria $X$ que en la población se distribuye con una distribución Normal(0.1,0.1). Y supongamos que estamos interesados en demostrar que la media de esa variable $X$ en la población es mayor a 0. Es decir:
Otras dos hipótesis alternativas son posibles; $H_a \ne 0$ ó $H_a < 0$. Podemos elegir cualquiera de éstas dependiendo de nuestro hipótesis de trabajo, determinando si el test será a una o dos colas. En este caso el test será a una cola y el p value para un valor observado $T$ será:
In [1]:
using Distributions
μ = 0.1
población = Normal(μ,0.1)
n = 5
muestra = rand(población, n)
Out[1]:
In [2]:
# H₀
μ₀ = 0
H₀ = TDist(n-1)
Out[2]:
In [3]:
x̄ = mean(muestra)
Out[3]:
In [4]:
s = std(muestra)
Out[4]:
In [5]:
T = (x̄-μ₀)/(s/sqrt(n))
Out[5]:
In [6]:
p = ccdf(H₀,T)
Out[6]:
Obtenemos el mismo resultado que utilizando OneSampleTTest de la biblioteca HypothesisTests de Julia:
In [7]:
using HypothesisTests
In [8]:
ttest = OneSampleTTest(muestra, 0)
Out[8]:
In [9]:
pvalue(ttest, tail=:right)
Out[9]:
El p value depende de si el test es a una o dos colas, en este caso, como la distribución es simétrica alrededor del 0, es obvio que el p value para el test a dos colas es el doble del p value que obtuvimos. Pero calcular el p value a dos colas no es trivial, principalmente puede ser un problema para distribuciones asimétricas (Kulinskaya, 2008; Agresti & Kateri, 2011). Pero en general se toma la idea de doblar el valor del P value menor al testear las dos opciones a una cola (es como si hiciéramos una corrección de Bonferroni al haber realizado dos pruebas).
Kulinskaya, Elena. "On two-sided p-values for non-symmetric distributions." arXiv preprint arXiv:0810.2124 (2008).
Agresti, Alan, and Maria Kateri. Categorical data analysis. Springer Berlin Heidelberg, 2011.
En particular, para una variable aleatoria $X$, que sigue una determinada distribución si $H_0$ es cierta, dado un estimador $x_0$, el p value es:
$P(X \geq x|H_0)$ prueba de una cola hacia la derecha (right-tail event)
In [10]:
p_right = pvalue(ttest, tail=:right)
Out[10]:
In [11]:
using Plots
pyplot(size=(400,200))
PDF = plot(t -> pdf(H₀,t), -10, T, legend=false)
plot!(PDF, t -> pdf(H₀,t), T, 10, fill=0)
CCDF = plot(t -> ccdf(H₀,t), -10, 10, legend=false)
scatter!(CCDF, [T], [p_right])
plot(PDF, CCDF)
Out[11]:
$P(X \leq x|H_0)$ prueba de una cola hacia la izquierda (left-tail event)
In [12]:
p_left = pvalue(ttest, tail=:left)
Out[12]:
In [13]:
PDF = plot(t -> pdf(H₀,t), T, 10, legend=false)
plot!(PDF, t -> pdf(H₀,t), -10, T, fill=0)
CDF = plot(t -> cdf(H₀,t), -10, 10, legend=false)
scatter!(CDF, [T], [p_left])
plot(PDF, CDF)
Out[13]:
$2 min{P(X \leq x|H_0), P(X \geq x|H_0)}$ para el test a dos colas (double-tailed event)
In [14]:
pvalue(ttest, tail=:both)
Out[14]:
In [15]:
2*min(p_left,p_right)
Out[15]:
In [16]:
PDF = plot(t -> pdf(H₀,t), -T, T, legend=false)
plot!(PDF, t -> pdf(H₀,t), -10, -T, fill=0)
plot!(PDF, t -> pdf(H₀,t), T, 10, fill=0)
CCDF = plot(t -> ccdf(H₀,t), -10, 10, legend=false)
scatter!(CCDF, [-T, T], [p_left,p_right])
plot(PDF, CCDF)
Out[16]:
Para facilitar ésto, la biblioteca HypothesisTests define la función pvalue para tomar directamente una distribución de probabilidad (no solamente un test estadístico), y calcular el p value para un determinado estadístico.
In [17]:
pvalue(TDist(n-1), T, tail=:right)
Out[17]:
El intervalo de confianza es un tipo de estimador por intervalo para un parámetro de la población. Se calcula a partir de una o más muestras de la población, y posee un nivel de confianza que determina qué tan frecuentemente una serie de intervalos estimados de la misma manera contiene al parámetro poblacional. Es importante notar que el intervalo de confianza calculado a partir de una única muestra puede no contener al parámetro a estimar, pero aún así es una buena estimación intervalar del mismo. Así un intervalo con una confianza del 95%, significa que de un número N de intervalos calculados, 0.95N contendrán al parámetro a estimar. No debe interpretarse nunca como un intervalo con una probabilidad del 95% de contener el parámetro.
El intervalo para una media se calcula utilizando el valor crítico de $t$ si la desviación estándar proviene de la muestra o $Z$ si se conoce la desviación estándar de la población. Los valores críticos son aquellos valores de $t$ o $Z$ que dejan un área de $\frac{\alpha}{2}$ hacia la izquierda del menor valor crítico, y el mismo hacia la derecha (dejando un área $1-\alpha$ en el intervalo).
Así el valor crítico, utilizando la distribución T de Student con $n-1$ grados de libertad, para una confianza del 95% ($\alpha = 5%$) puede calcularse con la función quantile de Distributions
In [18]:
t_left = quantile(TDist(n-1),0.05/2)
Out[18]:
In [19]:
t_right = quantile(TDist(n-1),1-0.05/2)
Out[19]:
Los valores son iguales en módulo dada la simetría de la distribución T de Student. Siendo $t$ ese módulo, el intervalo de confianza para la media de la población sería:
In [20]:
t = abs(t_left)
se = s/sqrt(n) # Error estándar de la media
Out[20]:
In [21]:
( x̄ - t*se, x̄ + t*se )
Out[21]:
Para facilitar ésto, la biblioteca HypothesisTests define la función ci.
In [22]:
confint(ttest) # usando ci de HypothesisTests sobre el T test
Out[22]:
In [23]:
n_experimentos = 50
intervalos = map(1:n_experimentos) do N
muestra = rand(población, n)
t_test = OneSampleTTest(muestra)
intervalo = confint(t_test, 0.1) # α = 0.1
[intervalo[1], intervalo[2]] # retorno un vector [inicio, fin]
end
Out[23]:
In [24]:
intervalos = hcat(intervalos...)
Out[24]:
In [25]:
y = hcat(map(N -> [N,N],1:n_experimentos)...)
Out[25]:
In [30]:
plot(intervalos, y, size=(600,600), legend=false, linecolor=:darkgrey, ylim=(-1,51), grid=false)
vline!([μ])
Out[30]:
In [31]:
contiene_μ = [ intervalos[1,i] ≤ μ ≤ intervalos[2,i] for i in 1:n_experimentos ]
Out[31]:
In [32]:
mean(contiene_μ) # ~ 90% ( α = 0.1, confianza = 1 - α = 0.9 )
Out[32]:
**Bootstrap** nos permite obtener estimaciones puntuales o por intervalos, en particular, intervalos de confianza, para parámetros poblacionales sin necesidad de conocer la distribución subyacente.
Básicamente la técnica de bootstrap se basa en hacer un muestreo aleatorio con reemplazo de nuestra muestra (remuestreo), calculando el parámetro de interés en cada muestra. La muestra de la que se parte debe ser una muestra representativa, en particular debe complir con ser una muestra de valores independientes e idénticamente distribuidos (iid). Si bien la distribución de la variable de interés en la población no forma parte de los supuestos del método, una implementación naïve puede generar resultados incorrectos con distribuciones con colas pesadas. Bootstrap también puede presentar problemas con distribuciones discretas y poco densas.
In [33]:
# Se hacen muchos remuestreos
medias = map(1:10_000) do i
# Se toma una muestra al azar (con reemplazo) del mismo tamaño que la muestra
resampling = sample(muestra, replace=true, n)
# Se calcula el estadístico de interés
mean(resampling)
end
# Bootstrap distribution of the statistic
histogram(medias, legend=false)
Out[33]:
In [47]:
quantile(medias,[0.05,0.95]) # α = 0.05 + 0.95 = 0.1, confianza = 1 - α = 0.9
Out[47]:
In [48]:
median(medias)
Out[48]:
A pesar de que ésta implementación naïve es super sencilla, la biblioteca Bootstrap de Julia ofrece herramientas para facilitar el proceso de bootstraping, así como implementaciones más complejas para reducir el error en algunas situaciones. Herramientas similares son ofrecidas por la biblioteca boot de R.
In [49]:
using Bootstrap
In [50]:
bootstrapped_mean = bootstrap(muestra, mean, BasicSampling(10_000))
Out[50]:
In [51]:
Bootstrap.ci(bootstrapped_mean, BasicConfInt(0.9)) # α = 0.1, confianza = 0.9
Out[51]:
Una ventaja de hacerlo con estas bibliotecas es que podemos tener intervalos de confianza con corrección del sesgo, BCa:
In [52]:
Bootstrap.ci(bootstrapped_mean, BCaConfInt(0.9)) # Bias-Corrected and accelerated (BCa)
Out[52]:
Muchas veces es necesario saber cuando los datos de una muestra se ajustan a una determinada distribución. En ese sentido uno se puede preguntar si la distribución observada (en general, para una muestra) es igual, o no, a una distribución esperada (en general, para la población). Para testear ésto, disponemos de las siguientes pruebas:
La prueba χ² se puede utilizar para testear la bondad de ajuste de cualquier distribución univariada para la cual se pueda calcular la CDF. Dado que requiere datos en grupos, por lo que puede utilizarse para testear distribución para datos categóricos (a diferencia de KS). Para poder utilizar el test en una distribución de probabilidad para datos continuos, debemos discretizar/agrupar los datos de la misma manera en que lo hacemos al calcular un histograma. Debemos tener en cuenta q el test es sensible a la elección de los bins. Una desventaja de este test es que requiere de un tamaño muestral de al menos 1000 para que la aproximación χ² sea válida.
In [53]:
población = Binomial(100, 0.055)
expected = Binomial(100, 0.05 ) # H₀: expected == observed
Out[53]:
In [54]:
observed = rand(población, 1000) # muestra
histogram(observed, legend=false, bins=-0.5:1:maximum(observed))
Out[54]:
In [55]:
using StatsBase
cdf_observed = ecdf(observed)
scatter( 0:1:maximum(observed), cdf_observed, alpha=0.7, label="observed")
scatter!(0:1:maximum(observed), x -> cdf(expected, x), alpha=0.7, label="expected")
Out[55]:
Dadas $k$ categorías, para la categoría $i$, $O_i$ representa la frecuencia observada para esa categoría, mientras $E_i$ representa la frecuencia esperada (nuestra hipótesis nula).
Ésta última, en el caso de haber generado bins para una variable continua, la calculamos usando la CDF evaluada en el extremo superior bin, menos su valor en el extremo inferior del bin. N es el tamaño de la muestra.
In [56]:
N = length(observed)
Oᵢ = Float64[ sum(observed .== x) for x in 0:1:100 ]
Eᵢ = Float64[ pdf(expected,x) for x in 0:1:100 ] .* N
Out[56]:
In [57]:
χ² = sum( ( (Oᵢ .- Eᵢ) .^ 2 ) ./ Eᵢ )
Out[57]:
In [58]:
k = 101 # Número de bins/categorías
c = 2 + 1 # Parámetros estimados para la distribución más uno
# k - c grados de libertad
# Se testea a una cola, a la derecha χ²
pvalue(Chisq(k-c), χ², tail=:right)
Out[58]:
In [59]:
H₀ = Chisq(k-c)
PDF = plot(x -> pdf(H₀,x), 0, χ², legend=false)
plot!(PDF, x -> pdf(H₀,x), χ², χ²+100, fill=0)
Out[59]:
In [61]:
ChisqTest(convert(Vector{Int}, Oᵢ), Eᵢ./sum(Eᵢ))
Out[61]:
La prueba de Kolmogorov-Smirnov es otra prueba de ajuste a una distribución que usa la función de densidad acumulada (CDF). En particular, su estadístico, $D$ mide la distancia entre la CDF de la distribución, que suponemos correcta en la hipótesis nula, y la ECDF de la muestra. En particular, el estadístico $D$ es el valor máximo de las diferencias. Si la función de densidad acumulada empírica es igual a la función de densidad de nuestra distribución a ajustar, el estadístico $D$ será 0 ($H_0$: $D = 0$).
Éste es un test paramétrico, por lo que hace pocas supuestos, pero debe tenerse en cuenta que la distribución que suponemos en nuestra $H_0$ no puede tener parámetros libres, debe estar completamente determinada y debe ser continua. Es sensible tanto a cambios en la ubicación como en la forma de la distribución, pero no tiene gran poder, por lo que puede requerir un tamaño muestral grande para rechazar correctamente la hipótesis nula.
In [62]:
población = Binomial(100, 0.05)
muestra = rand(población, 10)
ExactOneSampleKSTest(muestra, población)
Out[62]:
sleep dataset
Data which show the effect of two soporific drugs (increase in hours of sleep compared to control) on 10 patients.
- Extra: increase in hours of sleep
- Group: drug given
- ID: patient ID
In [63]:
using RDatasets
sleep_long = dataset("datasets","sleep")
Out[63]:
In [66]:
?DataFrames.unstack
Out[66]:
In [67]:
sleep_wide = unstack(sleep_long, :ID, :Group, :Extra)
Out[67]:
In [68]:
names!(sleep_wide, [:ID,:G1,:G2])
Out[68]:
In [70]:
using RCall
by(sleep_long, :Group) do df
extra = df[:Extra]
DataFrame(
P_AD = pvalue(OneSampleADTest(extra, Normal())),
P_Shapiro = rcopy(R"shapiro.test($extra)")[Symbol("p.value")]
)
end
Out[70]:
In [73]:
using StatPlots
boxplot(sleep_long, :Group, :Extra, legend=false)
Out[73]:
In [74]:
by(sleep_long, :Group) do df
var( df[:Extra] )
end
Out[74]:
Todos los test de igualdad de varianza tienen como hipótesis nula que las varianzas son iguales.
El test F de igualdad de varianzas tiene como supuesto que las distribuciones son normales.
In [75]:
R"var.test(Extra ~ Group, $sleep_long)"
Out[75]:
El test de Bartlett es sensible a las desviaciones de la normalidad, por lo que no es recomendable usarlo con distribuciones no normales.
In [76]:
R"bartlett.test(Extra ~ Group, $sleep_long)"
Out[76]:
El test de Levene es robusto ante distribuciones no normales.
In [79]:
# install.packages("car")
R"""
library(car)
leveneTest(Extra ~ Group, $sleep_long)
"""
Out[79]:
In [80]:
EqualVarianceTTest(sleep_wide[:G1], sleep_wide[:G2])
Out[80]:
In [81]:
UnequalVarianceTTest(sleep_wide[:G1], sleep_wide[:G2]) # Welch’s t-test
Out[81]:
In [82]:
OneSampleTTest(
convert(Vector{Float64}, sleep_wide[:G1]),
convert(Vector{Float64}, sleep_wide[:G2])
) # H₀: mean(x) - mean(y) = 0
Out[82]:
In [83]:
MannWhitneyUTest(sleep_wide[:G1], sleep_wide[:G2]) # H₀: P(x > y) = P(y > x)
Out[83]:
In [84]:
SignedRankTest(
convert(Vector{Float64}, sleep_wide[:G1]),
convert(Vector{Float64}, sleep_wide[:G2])
) # H₀: median(x - y) = 0
Out[84]:
In [85]:
ApproximateTwoSampleKSTest(sleep_wide[:G1], sleep_wide[:G2]) # H₀: x e y provienen de la misma distribución
Out[85]:
Then let W be the number of pairs for which yi − xi > 0. Assuming that H0 is true, then W follows a binomial distribution W ~ b(m, 0.5).
$m$ es el número de pares que muestran diferencias.
In [86]:
SignTest(sleep_wide[:G1], sleep_wide[:G2]) # H₀: P(x > y) = 0.5
Out[86]:
Es un caso especial de la prueba de chi-cuadrado. Esta prueba posee poco poder.
Se calcula la mediana para todos los datos ($x$ e $y$), estableciéndose dos grupos (por encima o por debajo de ésta). Luego se arma una tabla de contingencia con esta nueva variable versus las dos variables anteriores ($x$ e $y$).
In [87]:
median_xy = median(sleep_long[:Extra])
Out[87]:
In [88]:
sleep_long[:Side] = sleep_long[:Extra] .<= median_xy
Out[88]:
In [90]:
using FreqTables
In [91]:
tabla = freqtable(sleep_long, :Side, :Group)
Out[91]:
In [92]:
ChisqTest(tabla)
Out[92]:
Use Fisher's exact test when you have two nominal variables. You want to know whether the proportions for one variable are different among values of the other variable.
In [93]:
FisherExactTest(tabla...)
Out[93]:
Se utiliza para tablas de contingencia 2x2 de variables dicotómicas en medidas pareadas.
| X / Y | true | false |
|---|---|---|
| true | a | b |
| false | c | d |
Q sigue una distribución χ²con 1 grado de libertad para muestras grandes, es posible hacer una corrección por continuidad para conseguir una mejor aproximación a χ²:
In [125]:
X = rand(Binomial(1,0.6),100)
Y = rand(Binomial(1,0.4),100)
tabla = freqtable(X, Y)
Out[125]:
In [127]:
b = tabla[0,1]
Out[127]:
In [128]:
c = tabla[1,0]
Out[128]:
In [129]:
Q = ((b-c)^2)/(b+c)
Out[129]:
In [130]:
Qc = ((abs(b-c)-1)^2)/(b+c)
Out[130]:
In [131]:
pvalue(Chisq(1), Q, tail=:right)
Out[131]:
In [132]:
pvalue(Chisq(1), Qc, tail=:right)
Out[132]:
In [133]:
NamedArrays.setnames!(tabla, [1,2], 1)
NamedArrays.setnames!(tabla, [1,2], 2)
Out[133]:
In [134]:
R"mcnemar.test($tabla)"
Out[134]: