Diseños experimentales en R


1.1. Experimento básico

La dificultad de obtener datos experimentales en cualquier ciencia social dificulta observar efectos causales. Sin embargo, los diseños experimentales intentan aproximar la obtención de efectos causales a partir de datos observacionales, donde los factores muestran mayor variabilidad y menor control respecto a cuándo el experimento se realiza en el contexto de un laboratorio. En este marco, la variable resultado (o variable dependiente [y]) estaría siendo incidida por variables independientes que se pueden observar (adicional a nuestra variable de interés) y variables intervinientes no observadas, sobre la que debemos establecer supuestos para poder alcanzar un “efecto causal”.

1.2. Análisis estadístico preliminar

En la presente entrada se muestran algunos ejemplos de diseños experimentales utilizando R. De forma puntual se utiliza la data bwght del libro de Econometría de J. Wooldridge, disponible en el paquete wooldridge. El siguiente bloque de códigos adjunta la librería de ggplot2 para comparar el peso promedio de los niños al nacer (bwght) para dos grupos especifico, segmentados según hayan o no fumado las madres durante el embarazo (fuma). Note que es importante definir la variable categórica como un factor. La función is.factor nos permite verificar si una variable determinada está siendo considerada como un factor.

library(wooldridge)
library(ggplot2)

data(bwght)
attach(bwght)

bwght$fuma <- factor(packs>0, labels = c("No","Si"))

Es importante verificar el número de observaciones en cada categoría:

table(bwght$fuma)

  No   Si
1176  212

Un segundo antecedente es estudiar el peso promedio (o mediano) por grupo puede obtenerse en R de distintas maneras, aquí mostramos dos alternativas, la alternativa 1 está disponible en R (tapply) por default, la segunda utiliza el paquete dplyr de tidyverse. Observándose que el peso mediano es superior entre las mujeres que no fumaron durante el embarazo.

#Alt. 1
tapply(bwght$bwght, bwght$fuma, median)
No  Si
121 112

#Alt. 2
library(dplyr)
bwght %>%
  group_by(fuma) %>%
  summarise(median(bwght))

# A tibble: 2 x 2
  fuma  `median(bwght)`
  <fct>           <dbl>
1 No                121
2 Si                112

Además, una forma gráfica de obtener esta comparación es mediante el grafico de caja o boxplot, el mismo, adicional a indicarnos la mediana de la serie, nos permite ver la dispersión de la variable dependiente (bwght), lo que puede darnos una idea de la significancia en la diferencia de media. Otra relevancia del test es que permite verificar la presencia de valores atípicos.

ggplot(bwght, aes(y=bwght, x=fuma)) +
  geom_boxplot(outlier.colour="red") +
  theme_minimal()


Se puede observar que el peso promedio de los niños de madre que no fuman durante el embarazo es ligeramente superior al resto, pero con una amplia dispersión. En tal sentido, hay madres que fuman que tienen niños con un peso superior al de madres que no fuman. Esta comparación de promedios pudiese haber sido evidencia de un efecto causal en el marco de un experimento donde los tratamientos se asignan aleatoriamente a las unidades experimentales (Maguiña, 2018), pero como los datos no son experimentales, es necesario establecer el control de los factores observables (que estén disponibles en la base de datos) y establecer supuestos sobre aquellos factores no observados, resumida en aspectos como la equivalencia u homogeneidad entre los grupos. Es decir, no podemos a este nivel, indicar que fumar es la razón de la diferencia entre grupos.

1.1. Prueba t de comparación de media

Un test formal de comparación de media es el test t de comparación de media. Este test se realiza mediante el test t de comparación de media, que se realiza en R con la función t.test cuya hipótesis nula es que no existen diferencia entre la media de los grupos.

testt1 <- t.test(bwght$bwght~bwght$fuma,
                 alternative = "two.sided",
                 var.equal = FALSE,
                 conf.level = 0.95)
testt1

Welch Two Sample t-test

data:  bwght$bwght by bwght$fuma
t = 6.1743, df = 302.29, p-value = 2.138e-09
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  6.073641 11.756355
sample estimates:
mean in group No mean in group Si
        120.0612         111.1462

El p-valor obtenido claramente rechaza la hipótesis nula de igualdad en media. Recordemos que el p-valor representa el área de la función de densidad desde nuestro estadístico de prueba hasta el infinito, ahora, dado que la diferencia entre el estadístico crítico y el de prueba es que nos permite establecer una regla de rechazo de la nula, siempre rechazaremos la nula si el p valor es menor al nivel de significancia seleccionado.

Los paquetes moonBook y webr nos permiten obtener una versión grafica (posiblemente más fácil de interpretar) del test de hipótesis que hemos realizado anteriormente:

library(moonBook)
library(webr)
plot(testt1)

1.1. ¿Igualdad de la varianza?

Note que en el ejemplo anterior hemos supuesto varianzas distintas (var.equal = FALSE), basados en el boxplot. Una prueba formal que nos permite verificar si las varianzas son iguales es el test de Bartlett, cuya hipótesis nula es que las varianzas son iguales. En R, esta prueba se realiza mediante el comando bartlett.test:

library(car)
bartlett.test(bwght ~ nivelIng, data = bwght)

       Bartlett test of homogeneity of variances

data:  bwght by nivelIng
Bartlett's K-squared = 9.2027, df = 2, p-value = 0.01004

Adicionalmente, la prueba pwr pretende determinar el nivel de muestra necesario para alcanzar un poder determinado, o calcular el poder, dado un tamaño de muestra. Según los resultados, se necesitan muestra de al menos 83 individuos.

1.2. Test ANOVA

En el ejemplo, la mediana del peso de los niños cuyas madres han fumado durante el embarazo, es notablemente inferior al resto. Sin embargo, cuando necesitamos comparar más de dos grupos, como agregar el control de factores externos se utiliza el test ANOVA.

El próximo ejemplo realiza una comparación de media del peso al nacer entre tres grupos de ingresos. Para estos fines creamos la variable categórica nivelIng que crea tres niveles de ingresos:

bwght$nivelIng <- factor(cut(faminc, breaks = 3),
                         labels = c("bajo","medio","alto"))

ggplot(bwght, aes(y=bwght, x=nivelIng)) +
  geom_boxplot(outlier.colour="red") +

  theme_minimal()


También podemos usar las maravillas del tidyverse, para comparar el peso promedio de los niños al nacer para cada grupo de ingresos y observar que, aun controlando por nivel de ingresos, continúan observándose diferencias relativas en el peso de los niños según el consumo de cigarrillos de la madre. Ahora bien, note la diferencia relativa, representada con la variable diffrel, se reduce en la medida en aumentamos el nivel de ingresos del grupo, esto sería un indicador de los factores externos no observados, tales como la educación o el acceso atenciones médicas, estarían actuando en el problema.

bwght %>%
  group_by(nivelIng,fuma) %>%
  summarise(medianas = median(bwght)) %>%
  mutate(diffrel = (medianas/lag(medianas))-1)

# A tibble: 6 x 4
# Groups:   nivelIng [3]
  nivelIng fuma  medianas diffrel
  <fct>    <fct>    <dbl>   <dbl>
1 bajo     No        119  NA    
2 bajo     Si        110. -0.0714
3 medio    No        122  NA    
4 medio    Si        114. -0.0697
5 alto     No        124  NA    
6 alto     Si        120. -0.0282

Primero mostramos un ejemplo de comparación de media cuando tenemos varias categorías o clasificaciones, el test ANOVA (ANalysis Of VAriance), basados en independencia, normalidad y homocedasaticidad, que permite comparar distintos grupos (el test t visto anteriormente solo comparaba dos grupos). Siendo la hipótesis nula, que no existen diferencias estadísticas entre las medias de los grupos, el comando aov permite obtener la tabla Anova [referencia]:

anova_mod <- aov(bwght ~ nivelIng, data = bwght)
summary(anova_mod)

              Df Sum Sq Mean Sq F value   Pr(>F)   
nivelIng       2   8506    4253    10.4 3.27e-05 ***
Residuals   1385 566106     409                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Claramente se rechaza la hipótesis nula de no diferencias entre media.

1.1. Control de factores en el modelo ANOVA

Ahora, la mayor ventaja de esta metodología no es comparar varios grupos, sino poder introducir el control explícito de variables como el ingreso. El control de factores dentro de un modelo Anova se logra incluyendo estas variables de forma explícita en los modelos. Este análisis es importante:

a.       Cuando se necesita estudiar el efecto de varios factores.
b.       Para controlar el resto de variables relevantes y que han sido observadas. Recuerde que es necesario, como objetivo del diseño de la investigación, alcanzar que el residuo del modelo sea aleatorio, por tanto, independiente de x.  

Al realizar la comparación a lo interno de cada grupo de ingresos, estamos introduciendo el control explícito del nivel de ingreso (a menos que exista una alta dispersión a lo interno de cada grupo, de forma tal, que aun controlando por grupos de ingreso la desigualdad a lo interno de estos grupos pueda explicar diferencias importantes entro los individuos que lo componen).

El siguiente ejemplo muestra como agregar factores en R:

anova_mod1 <- aov(bwght ~ fuma + faminc, data = bwght)
summary(anova_mod1)

              Df Sum Sq Mean Sq F value   Pr(>F)   
fuma           1  14276   14276  35.527 3.19e-09 ***
faminc         1   3803    3803   9.464  0.00214 **
Residuals   1385 556533     402                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

1.2. Validar el modelo ANOVA

1.2.1.        Testear autocorrelación

La validez del modelo Anova depende fundamentalmente de que los residuos asuman un comportamiento aleatorio y normal. El Box-Pierce test testea independencia lineal (autocorrelación) entre las observaciones (una alternativa grafica se obtiene con los correlogramas acf(re)).

re <- anova_mod$residuals
Box.test(re)
Box-Pierce test

data:  re
X-squared = 2.58, df = 1, p-value = 0.1082

En el ejemplo anterior no es posible rechazar la hipótesis nula. Note, que en caso de que se rechace la hipótesis de no autocorrelación a un nivel del 5%, lo que indicaría cierta estructura de autocorrelación en los residuos. Esto podría ser resultado de alguna variable omitida, debido a que esta versión se corresponde con la del modelo de regresión simple. En términos llanos, esto indica que necesitaríamos intentar incluir nuevos modelos o modificar la forma funcional del modelo, para intentar no rechazar la nula.

1.2.2.        Testear normalidad

El gráfico qq-plot permite comparar la distribución empírica de nuestros datos, con aquella esperada en caso de que nuestra variable se distribuya como una normal. Usando el paquete tidyverse este se obtiene como:

ggplot(data.frame(re), aes(sample = re)) +
  stat_qq() +
  stat_qq_line() +
  theme_minimal()


Sin embargo, en muchos casos los test gráficos pueden no ser concluyentes. Para esto se usa el test shapiro.test, cuya hipótesis nula es que la variable de interés se comporta como una normal:

shapiro.test(re)

       Shapiro-Wilk normality test

data:  re
W = 0.97459, p-value = 6.276e-15

Cuando no se cumple con los supuestos de base, es recomendable usar un test no paramétrico.

kruskal.test(bwght$bwght~bwght$fuma)

       Kruskal-Wallis rank sum test

data:  bwght$bwght by bwght$fuma
Kruskal-Wallis chi-squared = 38.987, df = 1, p-value = 4.267e-10

1.1. Test de Tukey

Hasta aquí hemos usado un test paramétrico ANOVA para comparar media entre grupos, pero este no aporta información sobre las diferencias especificas entre los grupos. El test crea intervalos de confianzas y compara las distintas clases, para contrastar diferencias en media entre estas. De esta forma, podemos identificar a cuáles grupos específicos responden las posibles diferencias en media identificadas por el test Anova.

testtukey <- TukeyHSD(anova_mod)
testtukey

Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = bwght ~ nivelIng, data = bwght)

$nivelIng
               diff        lwr       upr     p adj
medio-bajo 3.912042  1.0809266  6.743158 0.0034789
alto-bajo  6.551258  2.9552061 10.147309 0.0000608
alto-medio 2.639215 -0.8809485  6.159379 0.1839291
 
El test anterior indica que no existen diferencia en el peso medio en los bebes de los grupos alto-medio, contrario a los demás grupos comparados, donde si se observan diferencias. Podemos graficar estos resultados utilizando la función plot. Esta prueba gráfica utiliza la lógica de un intervalo de confianza, si el cero está contenido en el intervalo de confianza, no podemos rechazar la nula de que esta diferencia sea igual a cero.
  
plot(testtukey)


1.1. Anova factor

El factor permite estimar diversas combinaciones de un grupo de variables. La forma funcional de esta especificación incluye interacciones entre estás.

anova_mod2 <- aov(bwght ~ fuma * male * white, data = bwght)
summary(anova_mod2)

                Df Sum Sq Mean Sq F value   Pr(>F)   
fuma               1  14276   14276  36.015 2.50e-09 ***
male               1   2481    2481   6.259   0.0125 * 
white              1   9333    9333  23.545 1.36e-06 ***
fuma:male          1    544     544   1.372   0.2417   
fuma:white         1    507     507   1.280   0.2582   
male:white         1     36      36   0.090   0.7644   
fuma:male:white    1    425     425   1.072   0.3006   
Residuals       1380 547011     396                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

En este caso, los estadísticos de prueba son distintos. Note que una vez controlada la raza y el sexo del niño, no se observan diferencias significativas en la media del peso de los niños. Ahora bien, los estadísticos de prueba ahora son distintos, estos se obtienen mediante:

drop1(anova_mod2, ~., test="F")
Single term deletions

Model:
bwght ~ fuma * male * white
                Df Sum of Sq    RSS    AIC F value    Pr(>F)   
<none>                       547011 8311.5                     
fuma             1      16.3 547027 8309.6  0.0412 0.8392641   
male             1    1479.6 548490 8313.3  3.7327 0.0535627 . 
white            1    5763.7 552774 8324.1 14.5406 0.0001432 ***
fuma:male        1     863.4 547874 8311.7  2.1781 0.1402158   
fuma:white       1     941.5 547952 8311.9  2.3752 0.1235033   
male:white       1     186.2 547197 8310.0  0.4699 0.4931618   
fuma:male:white  1     425.0 547436 8310.6  1.0722 0.3006339   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Referencias

-          Calvo, J.; Polazón, F.; Carreño, P. (2016). Diseño experimental y análisis estadístico. Universidad de Murcia. Consultado 27/2/2020. Disponible en: https://aulavirtual.um.es/access/content/group/COLLAB-x5vsnmqzqac8cs9edm6rhxm/DISE%C3%91O%20EXPERIMENTAL%20Y%20AN%C3%81LISIS%20DE%20DATOS/DEyAE%20guion.pdf

-          Maguiña, M. (2018). Introducción a los diseños experimentales. Consultado 27/2/2020. Disponible en: https://rpubs.com/misael892/DCAyBCA.

-          Mendiburu, F. (2007). Diseño y Análisis de experimentos Aplicados en investigación Agrícola


Entradas populares de este blog

Valores perdidos (NA) en R: identificación y tratamiento (I)

Ejemplos de tablas en Stata

Métodos y técnicas de Investigación Económica