[1] 101.3561[1] 13.69224Estadística aplicada con R
Módulo 4: Herramientas Estadísticas
Muestreo simple sin reemplazo
Tomado de Learning Statistics with R
Muestreo simple con reemplazo
Tomado de Learning Statistics with R
Distribución de muestreo
La idea de poder medir enormes números de individuos es irreal.
Sin embargo, consideremos el siguiente escenario: si en lugar de medir el IQ de 100000 personas, repito el experimento una y otra vez pero en grupos de 5, puedo observar que la distribución de las medias aritméticas de todos estos experimentos adopta la forma de una distribución normal
Distribución de muestreo
Esta distribución toma el nombre distribución de muestreo de la media aritmética.
Lo que nos demuestra es que, incluso ante reducidos números de observaciones en una muestra, la media aritmética de esta muestra (\(\overline{X}\)) estará próxima a la media aritmética verdadera de la población (\(\mu\)).
Teorema del límite central
El teorema del límite central establece que siempre que el número de observaciones sea lo suficientemente grande, la distribución de muestreo de la media aritmética tenderá a ser normal independientemente de si la distribución de las observaciones es normal o no.
Ejemplo: supongamos que el ancho del caparazón de una especie de tortugas está comprendido entre 4 y 10 centímetros. En otras palabras, si observamos al azar una tortuga de esta población, sabemos que el ancho del caparazón estará en este rango. En términos de una distribución, podemos decir que el ancho del caparazón en una población de 1000 tortugas podría verse de esta manera:
Teorema del límite central
Pruebas de dos colas
\[\begin{align} H_0& : \mu_c- \mu_r = \theta \\ H_a& : \mu_c- \mu_r \neq \theta \end{align}\]
Imagen tomada de UCLA: Advanced Research Computing
Pruebas de una cola
\[\begin{align} H_0& : \mu_c- \mu_r \ge \theta \\ H_a& : \mu_c- \mu_r < \theta \end{align}\]
Imagen tomada de UCLA: Advanced Research Computing
Pruebas de una cola
\[\begin{align} H_0& : \mu_c- \mu_r \le \theta \\ H_a& : \mu_c- \mu_r > \theta \end{align}\]
Imagen tomada de UCLA: Advanced Research Computing
Antes de continuar
Gráfico Q-Q
Supuestos en la práctica
Usemos las pruebas con datos reales, esta vez con la tabla de datos
crimede la libreríaUsing R.Esta tabla de datos contiene registros de las tasas de crimen (# de reportes/100000 habitantes) en 50 estados en los Estados Unidos correspondiente a los años 1983 y 1993.
Sin enfocarnos por el momento en que prueba t específica usaremos, limitémonos a chequear la normalidad y la homogeneidad de varianzas entre las tasas de crimen registradas en 1983 y 1993.
Homogeneidad de las varianzas
F test to compare two variances
data: crime$y1983 and crime$y1993
F = 0.49163, num df = 50, denom df = 50, p-value = 0.01342
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2806163 0.8613069
sample estimates:
ratio of variances
0.4916266 
Transformación de variables
Raíz cuadrada
Logaritmo
Homogeneidad de las varianzas con transformación logarítmica
F test to compare two variances
data: log(crime$y1983) and log(crime$y1993)
F = 0.84048, num df = 50, denom df = 50, p-value = 0.5412
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.4797393 1.4724831
sample estimates:
ratio of variances
0.8404808 
Tabla de contingencia \(n \times n\)
Dependencia entre las filas y las columnas
Tabla de contingencia \(2 \times 2\)
La variable B (o respuesta) contendrá como posibles resultados “éxito” o “fracaso”.
La variable A (o explanatoria) posee las clases que identifican los grupos cuya probabilidad es sujeto de comparación.
El dataset de recursos por depredación
Gráficos de barras
Gráficos de barras
p1 <- ggbarplot(data = ranas,
x = "Pred",
y = "Age.FromEmergence",
add = "mean_se",
fill = "Pred")
p1 + labs(title = "Gráfico de barras de edad desde oviposición",
subtitle = "Datos del estudio de Prof. Touchon",
caption = "Gráfica propia",
x = "Predadores",
y = "Edad desde oviposición",
color = "Predador")
Gráficos de caja y bigote
Gráficos de densidad
Matriz de dispersión
Diagnósticos del ANOVA
Diagnósticos del ANOVA
Residuos vs. Valores ajustados
En este plot podemos evidenciar departuras del supuesto de la homocedasticidad. Idealmente, la línea roja que se muestra debería ir a lo largo de la horizontal en la coordenada cero del eje y (sobre la línea entrecortada).
Diagnósticos del ANOVA
Gráfico Q-Q
A diferencia del gráfico Q-Q que vimos para las pruebas t, en el eje y de este mismo gráfico para el ANOVA (y regresión lineal) se representan los residuos estandarizados. La interpretación es la misma: idealmente los puntos deberían ir a lo largo de la diagonal. Cuando no es así, evidencia una violación del supuesto de la normalidad.
Diagnósticos del ANOVA
Raíz cuadrada de los residuos estandarizados vs. Valores ajustados
Similar al primer diagnóstico, en el caso del ANOVA, nos da una idea de posibles departuras de la homogeneidad de las varianzas. La línea roja idealmente debería ser completamente recta.
Diagnósticos del ANOVA
Residuos vs. Apalancamiento
Aquellos puntos que estén etiquetados con números son mostrados como posibles outliers bajo dos criterios:
Están por fuera de los límites de la regla del rango intercuartílico (IQR), y
Marcados como outliers con influencia de apalancamiento mediante la prueba de Cook (distancia de Cook).
El segundo criterio es un argumento sólido para remover outliers.
Transformaciones
¿Qué tan distintos?
Preguntas naturales seguidas de estos resultados son: ¿qué tan distintos son los tratamientos entre sí?, ¿puedo acaso ordernarlos de mayor a menor?, ¿existen pares de tratamientos que son iguales?
Podemos comenzar con una ayuda visual (ligeramente distinto a nuestro AED, usando la transformación logarítmica)
Antes de continuar
Gráficos de comparaciones múltiples
ranas$Pred <- factor(ranas$Pred, levels = c("C", "NL", "L"))
bxplot <- ggboxplot(ranas, x = "Pred",
y = "Age.FromEmergence",
color = "Pred")
hsdvals <- emmeans_test(log.Age.FromEmergence ~ Pred,
data = ranas,
p.adjust.method = "mvt")
hsdvals <- add_significance(hsdvals,
p.col = "p.adj",
output.col = "p.adj.signif",
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", "ns"))
hsdvals <- hsdvals %>% add_xy_position(x = "Pred")
bxplot + stat_pvalue_manual(hsdvals, y.position = c(120, 130, 140))
Gráficos de comparaciones múltiples
library(stringr)
gruposvals <- as.data.frame(cld(ph1))
gruposvals$Pred <- factor(gruposvals$Pred, levels = c("C", "NL", "L"))
ggplot(gruposvals,
aes(x = Pred,
y = response,
fill = Pred)) +
geom_bar(stat = "identity",
show.legend = F,
color = "black")+
geom_errorbar(aes(ymin = response - SE,
ymax = response + SE),
width=0.2)+
geom_text(aes(label=str_trim(.group),
y = response+SE, vjust=-0.5))
ANOVA de múltiples vías en R
ANOVA de múltiples vías en R. Opc. 2
El inconveniente de usar la corrección de Welch es que esta se puede aplicar solo a ANOVAs de una vía.
Entonces, de manera similar a lo hicimos para poder analizar por Kruskall-Wallis un diseño factorial, tenemos que crear una variable dummy combinando los dos factores para correr un ANOVA de una vía.
Debemos tener en cuenta que al realizar esto, reducimos poder de la prueba estadística.
ranas$log.Age.DPO <- log(ranas$Age.DPO)
ranas$tratamiento <- paste(ranas$Pred, ranas$Res, sep = "-")
ranas$tratamiento <- factor(ranas$tratamiento, levels = c("C-Lo", "NL-Lo", "C-Hi", "L-Lo", "NL-Hi", "L-Hi"))
anova_op2 <- ranas %>%
welch_anova_test(log.Age.DPO ~ tratamiento)
anova_op2# A tibble: 1 × 7
.y. n statistic DFn DFd p method
* <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
1 log.Age.DPO 78 17.2 5 24.0 0.000000305 Welch ANOVA- Ahora podemos realizar las comparaciones de HSD Tukey (con corrección de Games-Howell)
games_comp <- ranas %>% games_howell_test(log.Age.DPO ~ tratamiento)
games_comp <- add_significance(games_comp,
p.col = "p.adj",
output.col = "p.adj.signif",
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", "ns"))
games_comp# A tibble: 15 × 8
.y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 log.Age.DPO C-Lo NL-Lo -0.155 -0.772 0.462 0.936 ns
2 log.Age.DPO C-Lo C-Hi -0.440 -0.700 -0.180 0.000261 ***
3 log.Age.DPO C-Lo L-Lo -0.485 -0.754 -0.217 0.0000932 ***
4 log.Age.DPO C-Lo NL-Hi -0.493 -0.821 -0.166 0.002 **
5 log.Age.DPO C-Lo L-Hi -0.643 -0.877 -0.408 0.00000114 ***
6 log.Age.DPO NL-Lo C-Hi -0.285 -0.902 0.331 0.545 ns
7 log.Age.DPO NL-Lo L-Lo -0.331 -0.947 0.286 0.419 ns
8 log.Age.DPO NL-Lo NL-Hi -0.339 -0.961 0.284 0.433 ns
9 log.Age.DPO NL-Lo L-Hi -0.488 -1.11 0.133 0.129 ns
10 log.Age.DPO C-Hi L-Lo -0.0452 -0.258 0.168 0.986 ns
11 log.Age.DPO C-Hi NL-Hi -0.0531 -0.351 0.244 0.988 ns
12 log.Age.DPO C-Hi L-Hi -0.202 -0.362 -0.0425 0.008 **
13 log.Age.DPO L-Lo NL-Hi -0.00797 -0.310 0.294 1 ns
14 log.Age.DPO L-Lo L-Hi -0.157 -0.334 0.0196 0.099 .
15 log.Age.DPO NL-Hi L-Hi -0.149 -0.440 0.142 0.448 ns 
Gráficos de interacción en R
Diagnósticos de la regresión lineal
Diagnósticos de la regresión lineal
Residuos vs. Valores ajustados
En el ANOVA vimos como este plot sugería departuras de la homocedasticidad. En el caso de la regresión lineal, los residuos al no estar agrupados en categorías presentan mayor información sobre este supuesto. Adicionalmente, curvaturas en la línea roja evidencian también departuras de la linearidad. Esto quiere decir que la relación entre las variables no es completamente lineal. A veces esto puede corregirse con transformaciones.
Diagnósticos de la regresión lineal
- Es ideal. (b) Es indicativo de no linearidad. (c) Evidencia de heterodasticidad. (d) Evidencia de una tendencia temporal
Diagnósticos de la regresión lineal
Gráfico Q-Q
Misma explicación que antes
Diagnósticos de la regresión lineal
Raíz cuadrada de los residuos estandarizados vs. Valores ajustados
Misma explicación que antes
Diagnósticos de la regresión lineal
Residuos vs. Apalancamiento
Aquellos puntos que estén etiquetados con números son mostrados como posibles outliers bajo dos criterios:
Están por fuera de los límites de la regla del rango intercuartílico (IQR), y
Marcados como outliers con influencia de apalancamiento mediante la prueba de Cook (distancia de Cook).
El segundo criterio es un argumento sólido para remover outliers.
Transformación de datos
Al igual que en el ANOVA, podemos recurrir a transformaciones de datos que nos permitan aliviar varios de los supuestos no cumplidos de la regresión lineal.
En el presente caso, esto lo logramos al utilizar el logaritmo natural de las dos variables del modelo
Call:
lm(formula = log(SVL.final) ~ log(Age.DPO), data = ranas)
Residuals:
Min 1Q Median 3Q Max
-0.182298 -0.049070 0.000516 0.038342 0.159057
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.44001 0.09211 37.345 < 2e-16 ***
log(Age.DPO) -0.11624 0.02232 -5.208 1.58e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.06244 on 76 degrees of freedom
Multiple R-squared: 0.263, Adjusted R-squared: 0.2533
F-statistic: 27.12 on 1 and 76 DF, p-value: 1.581e-06Interpretación de la regresión lineal
\[ y = mx + b \]
Call:
lm(formula = log(SVL.final) ~ log(Age.DPO), data = ranas)
Residuals:
Min 1Q Median 3Q Max
-0.182298 -0.049070 0.000516 0.038342 0.159057
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.44001 0.09211 37.345 < 2e-16 ***
log(Age.DPO) -0.11624 0.02232 -5.208 1.58e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.06244 on 76 degrees of freedom
Multiple R-squared: 0.263, Adjusted R-squared: 0.2533
F-statistic: 27.12 on 1 and 76 DF, p-value: 1.581e-06- Por cada incremento en una unidad del logaritmo de
Age.DPO, tenemos una unidad en descenso del logaritmo deSVL.final
Interpretación de la regresión lineal
Sin embargo una interpretación en la escala logarítmica no es completamente entendible, al menos para alguien ajeno al análisis que realizamos.
Podríamos hacer la retransformación de las variables a sus unidades originales y generar un gráfico de las predicciones. A la final se reduce a matemática básica.
Afortunadamente la librería
ggeffectsnos salva de ese dilema!
$Age.DPO
# Predicted values of SVL.final
Age.DPO | Predicted | 95% CI
------------------------------------
35 | 20.63 | [20.05, 21.23]
50 | 19.79 | [19.46, 20.13]
60 | 19.38 | [19.11, 19.65]
75 | 18.88 | [18.57, 19.20]
90 | 18.48 | [18.08, 18.90]
105 | 18.16 | [17.66, 18.67]
120 | 17.88 | [17.30, 18.48]
145 | 17.49 | [16.79, 18.22]
attr(,"class")
[1] "ggalleffects" "list"
attr(,"model.name")
[1] "lm6"
El mágico \(R^2\)
¡Gracias!
