Prueba de hipótesis

Realiza los siguientes ejercicios y genera un documento en RMarkdown con el código, los resultados y comentarios con tu interpretación. Recuerda utilizar análisis exploratorios y de supuestos estadísticos para determinar qué tipo de prueba de hipótesis es adecuada en cada caso. Justifica tu elección e interpreta los resultados.

Ejercicio 1

¿Reducir los servicios o aumentar los impuestos? En estos días, ya sea a nivel local, estatal o nacional, el gobierno a menudo enfrenta el problema de no tener suficiente dinero para pagar los diversos servicios que brinda. Una forma de abordar este problema es aumentar los impuestos. Otra forma es reducir los servicios. ¿Cual preferirías? Cuando la Encuesta de Florida preguntó recientemente a una muestra aleatoria de 1200 floridanos, el 52% (624 de los 1200) dijo que aumentaría los impuestos y el 48% dijo que reduciría los servicios. Determina si quienes están a favor de aumentar los impuestos en lugar de reducir los servicios son mayoría o minoría de la población.

datos <- data.frame(name=c("impuestos","servicios"), count=c(624, 1200-624))

respuestas <- rep(datos$name, times = datos$count)

df_respuestas <- data.frame(respuesta = respuestas)

str(df_respuestas)

## 'data.frame':    1200 obs. of  1 variable:
##  $ respuesta: chr  "impuestos" "impuestos" "impuestos" "impuestos" ...

library(gmodels)

CrossTable(df_respuestas$respuesta)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  1200 
## 
##  
##           | impuestos | servicios | 
##           |-----------|-----------|
##           |       624 |       576 | 
##           |     0.520 |     0.480 | 
##           |-----------|-----------|
## 
## 
## 
## 

Al explorar los datos se determina que se trata de una sola muestra con datos categóricos. Como el objetivo es comparar proporciones de la misma muestra, la prueba de hipótesis a utilizar es la prueba Chi-cuadrado bilateral.

library(ggstatsplot)

ggpiestats(df_respuestas, x=respuesta)

chisq.test(table(df_respuestas$respuesta), p = c(1/2, 1/2))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(df_respuestas$respuesta)
## X-squared = 1.92, df = 1, p-value = 0.1659

Cuando se obtiene el resultado de la prueba Chi-cuadrado bilateral, el p-valor no es inferior al 0.05 por lo que no podemos rechazar la hipótesis nula de igualdad, es decir que no se encontró evidencia de que quienes están a favor de aumentar los impuestos en lugar de reducir los servicios son mayoría o minoría de la población.

También, al ver el tamaño del efecto según la C de Pearson (0.04), es un tamaño del efecto muy pequeño, este resultado se obtiene en un intervalo de confianza entre 0.00 y 0.10 a largo plazo, en aproximadamente 95 de 100 veces que se repita el experimento.

Ejercicio 2

Se quiere evaluar un estudio de gemelos del mismo sexo donde un gemelo había tenido una condena penal. Se recopiló la siguiente información: si el hermano también había tenido una condena penal y si los gemelos eran gemelos monocigóticos (idénticos) o dicigóticos (no idénticos). Los estudios de gemelos como este se han utilizado a menudo para investigar los efectos de la “naturaleza versus crianza”. Queremos contrastar si la proporción de condenados es mayor para los gemelos monocigóticos que para los dicigóticos.

datos2 <- data.frame(
  Tipo = factor(c(rep("Dizygotic", 2 + 15), 
                  rep("Monozygotic", 10 + 3))),   
  Status = factor(c(rep("Convicted", 2),          
                    rep("Not_convicted", 15),     
                    rep("Convicted", 10),         
                    rep("Not_convicted", 3)))      
)

datos2

##           Tipo        Status
## 1    Dizygotic     Convicted
## 2    Dizygotic     Convicted
## 3    Dizygotic Not_convicted
## 4    Dizygotic Not_convicted
## 5    Dizygotic Not_convicted
## 6    Dizygotic Not_convicted
## 7    Dizygotic Not_convicted
## 8    Dizygotic Not_convicted
## 9    Dizygotic Not_convicted
## 10   Dizygotic Not_convicted
## 11   Dizygotic Not_convicted
## 12   Dizygotic Not_convicted
## 13   Dizygotic Not_convicted
## 14   Dizygotic Not_convicted
## 15   Dizygotic Not_convicted
## 16   Dizygotic Not_convicted
## 17   Dizygotic Not_convicted
## 18 Monozygotic     Convicted
## 19 Monozygotic     Convicted
## 20 Monozygotic     Convicted
## 21 Monozygotic     Convicted
## 22 Monozygotic     Convicted
## 23 Monozygotic     Convicted
## 24 Monozygotic     Convicted
## 25 Monozygotic     Convicted
## 26 Monozygotic     Convicted
## 27 Monozygotic     Convicted
## 28 Monozygotic Not_convicted
## 29 Monozygotic Not_convicted
## 30 Monozygotic Not_convicted

Luego de explorar los datos y ver los objetivos del análisis, se trata de dos grupos de entidades independientes de los cuales hay que comparar las proporciones de dos de sus características descritas por dos variables categóricas nominales. La prueba de hipótesis a utilizar es para comparar proporciones de dos muestras independientes ya que los datos vienen de entidades que no están relacionadas. Las pruebas de Chi-cuadrado para dos muestras independientes y la prueba exacta de Fisher en este caso son la elección.

CrossTable(datos2$Tipo, datos2$Status, fisher = TRUE, prop.r = FALSE, prop.c = TRUE, prop.t = FALSE, prop.chisq = FALSE, expected = TRUE)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |              Expected N |
## |           N / Col Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  30 
## 
##  
##              | datos2$Status 
##  datos2$Tipo |     Convicted | Not_convicted |     Row Total | 
## -------------|---------------|---------------|---------------|
##    Dizygotic |             2 |            15 |            17 | 
##              |         6.800 |        10.200 |               | 
##              |         0.167 |         0.833 |               | 
## -------------|---------------|---------------|---------------|
##  Monozygotic |            10 |             3 |            13 | 
##              |         5.200 |         7.800 |               | 
##              |         0.833 |         0.167 |               | 
## -------------|---------------|---------------|---------------|
## Column Total |            12 |            18 |            30 | 
##              |         0.400 |         0.600 |               | 
## -------------|---------------|---------------|---------------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  13.03167     d.f. =  1     p =  0.0003062666 
## 
## Pearson's Chi-squared test with Yates' continuity correction 
## ------------------------------------------------------------
## Chi^2 =  10.45814     d.f. =  1     p =  0.001221099 
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio:  0.04693661 
## 
## Alternative hypothesis: true odds ratio is not equal to 1
## p =  0.0005367241 
## 95% confidence interval:  0.003325764 0.3631823 
## 
## Alternative hypothesis: true odds ratio is less than 1
## p =  0.0004651809 
## 95% confidence interval:  0 0.2849601 
## 
## Alternative hypothesis: true odds ratio is greater than 1
## p =  0.9999845 
## 95% confidence interval:  0.00501032 Inf 
## 
## 
## 

ggbarstats(data = datos2, x = Tipo, y = Status, bf.message = FALSE)

Al ajustar la prueba Chi-cuadrado vemos que se rechaza la hipótesis nula de igualdad entre los grupos, ya que el p-valor < 0.05, tanto con el método tradicional como con la corrección de continuidad de Yates. Tambien, cuando vemos el resultado de la prueba exacta de Fisher se obtuvo un p-valor inferior. Esto sugiere que hay una relación entre el tipo de gemelos y el estatus criminal, donde en el grupo de convictos hay una mayor proporción de gemelos monocigóticos.

La puntuación obtenida para el tamaño del efecto fue de 0.64 según la V de Cramer, para un intervalo de confianza entre 0.24 y 1.0, lo que corresponde a un tamaño del efecto grande. En ese mismo sentido, en la prueba de Fisher se obtuvieron probabilidades relativas (odds ratio) de 0.047, menores que 1, a favor de la hipótesis de que existe una diferencia donde los condenados son mayor proporción en gemelos monocigóticos. Por lo que podemos concluir que existe una asociación entre ambas variables y se trata de una asociación fuerte.

Ejercicio 3

Vamos a evaluar si existe una relación entre el nivel educativo y el número abortos inducidos. La base de datos infert corresponde a un estudio de caso-control donde la variable “Education” está formada por 3 categorías (0 = 0-5 años, 1 = 6-11 años, 2 = 12+ años); y la variable “number of prior induced abortions” también (0 = 0, 1 = 1, 2 = 2 o más abortos inducidos).

Para acceder a los datos escribe en la consola de R:

data(infert, package = "datasets")

tab <- table(infert$education, infert$induced)

tab

##          
##            0  1  2
##   0-5yrs   4  2  6
##   6-11yrs 78 27 15
##   12+ yrs 61 39 16

Luego de ver y explorar lo datos utilizando una tabla de contigencias, se trata de dos variables categóricas ordinales donde se quiere determinar si existe o no alguna relación entre el nivel educativo y el número de abortos inducidos, en grupos que no están relacionados, por lo que la prueba de hipótesis a utilizar es para comparar proporciones de más de dos muestras independientes.

En este caso la selección fue la prueba de Fisher para comparaciones múltiples, ya que las frecuencias son bajas y esta opción resulta mas sencilla y adecuada.

library(RVAideMemoire)

fisher.multcomp(tab)

## 
##         Pairwise comparisons using Fisher's exact test for count data
## 
## data:  tab
## 
##                    0:1     0:2     1:2
## 0-5yrs:6-11yrs  0.7319 0.04040 0.10024
## 0-5yrs:12+ yrs  1.0000 0.05306 0.05306
## 6-11yrs:12+ yrs 0.1002 0.70712 0.70712
## 
## P value adjustment method: fdr

Después de aplicar la prueba exacta de Fisher que utiliza el ajuste por comparaciones múltiples (FDR), solo se encontró una diferencia estadísticamente significativa entre las mujeres con 0–5 años de educación y aquellas con 6–11 años de educación, específicamente al comparar la frecuencia de 0 abortos inducidos frente a 2 o más abortos inducidos (p = 0.0404). Este resultado sugiere que existe una asociación entre las variables pero solo para estos grupos en particular.

Otra comparación que presento una diferencia estadísticamente significativa marginal fue al comparar el grupo de educación de 0-5 años y el de 12 y más en los renglones 0 y 2 o más abortos (p = 0.053). Estos grupos educacionales también presentaron un p-valor semejante (p = 0.053) al comparar las categorías 1 y 2 o más abortos.

library(dplyr)

infert2 <- infert %>%
     filter(education %in% c("0-5yrs", "6-11yrs"),
         induced %in% c(0, 2)) %>% 
  select(education, induced)

rstatix::cramer_v(infert2$education, infert2$induced)

## [1] 0.2817095

Según la V de Cramer, el tamaño del efecto al comparar las diferencias entre los grupos de 0-5 años de educación y 6 a 11 años, en las categorías 0 abortos frente a 2 o más, se trata de un tamaño del efecto pequeño (0.28).

infert3 <- infert %>%
  filter(education %in% c("0-5yrs", "12+ yrs"),
    induced %in% c(0, 2)) %>%
      select(education, induced)
    
rstatix::cramer_v(infert3$education, infert3$induced)

## [1] 0.2463469

infert4 <- infert %>%
  filter(
    education %in% c("0-5yrs", "12+ yrs"),
    induced %in% c(1, 2)
  ) %>%
      select(education, induced)

rstatix::cramer_v(infert4$education, infert4$induced)

## [1] 0.2706417

El tamaño del efecto también fue pequeño para las comparaciones que obtuvieron diferencias significativas marginales.

Ejercicio 4

Utiliza los datos “Arthritis”, del paquete “vcd”, sobre un ensayo clínico de doble ciego que investiga un nuevo tratamiento para la artritis reumatoide.Tenemos información de 84 observaciones de 5 variables: la identificación del paciente (ID), el tratamiento (Treatment: Placebo, Treated), el sexo (Sex: Female, Male), la edad (Age) y la mejoría (Improved: None, Some, Marked). Para el grupo tratamiento, queremos comparar las edades de los pacientes que no mostraron mejoría con los que sí tuvieron una marcada mejoría.

library(vcd)

datos4 <- Arthritis %>%
   filter(Treatment == "Treated" & Improved %in% c("None", "Marked")) %>%
   select(Age, Improved) %>%
   group_by(Improved)

head(datos4)

## # A tibble: 6 × 2
## # Groups:   Improved [2]
##     Age Improved
##   <int> <ord>   
## 1    29 None    
## 2    30 None    
## 3    32 Marked  
## 4    46 Marked  
## 5    58 Marked  
## 6    59 None

El análisis a realizar es sobre las edades de dos grupos independientes, por lo que la prueba a utilizar en este caso es para datos numéricos.

library(rstatix)

datos4 %>%
     group_by(Improved) %>%
     get_summary_stats(Age, type = "mean_sd")

## # A tibble: 2 × 5
##   Improved variable     n  mean    sd
##   <ord>    <fct>    <dbl> <dbl> <dbl>
## 1 None     Age         13  49.8 16.8 
## 2 Marked   Age         21  56.8  9.02

library(ggpubr)

ggboxplot(x = "Improved", y = "Age",
           data = datos4,
           add = c("mean"),
           add.params = list(color = 'red'))

Luego de explorar los datos, la media del grupo Marked es superior a la del grupo None. Las desviaciones estandar reflejan una mayor diferencia, donde el grupo None presenta casi el doble de variabilidad en los datos.

datos4 <- as.data.frame(datos4)

datos4 %>%
     group_by(Improved) %>%
     identify_outliers(Age)

## # A tibble: 2 × 4
##   Improved   Age is.outlier is.extreme
##   <ord>    <int> <lgl>      <lgl>     
## 1 Marked      32 TRUE       TRUE      
## 2 Marked      41 TRUE       FALSE

Al evaluar supuestos generales para saber cual sería la prueba mas apropiada en este caso, vemos que tenemos valores atípicos y extremos, por lo que descartamos la prueba T para 2 muestras independientes.

datos4 %>%
     levene_test(Age ~ Improved)

## # A tibble: 1 × 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     1    32      6.68 0.0145

Continuando con la revisión de los supuestos, se rechaza la hipótesis de igualdad de varianza, aunque las pruebas restantes en este caso requieren revisar nuevamente este supuesto más adelante.

datos4 %>%
  group_by(Improved) %>%
  shapiro_test(Age)

## # A tibble: 2 × 4
##   Improved variable statistic      p
##   <ord>    <chr>        <dbl>  <dbl>
## 1 None     Age          0.877 0.0659
## 2 Marked   Age          0.903 0.0395

ggqqplot(datos4,
         x = "Age",
         facet.by = "Improved")

Al revisar la normalidad en la distribución de las edades en ambos grupos, según la prueba de normalidad de Shapiro-Wilk, se puede rechazar la hipótesis nula de no-normalidad en el grupo Marked, aunque el grupo None también presenta una nivel marginal de significación estadística.

En el gráfico quantil-quantil confirmamos que ninguna de las variables presenta una distribución normal.

Concluimos que los datos tampoco cumplen con el supuesto de normalidad. Procedemos entonces a revisar los supuestos para poder realizar la prueba robusta de Yuen para 2 muestras independientes.

datos4 %>%
    group_by(Improved) %>%
    filter(between(Age,
                   quantile(Age, 0.1),
                   quantile(Age, 0.9))) %>%
    get_summary_stats(Age, type = "mean_sd")

## # A tibble: 2 × 5
##   Improved variable     n  mean    sd
##   <ord>    <fct>    <dbl> <dbl> <dbl>
## 1 None     Age          9  50.9 13.6 
## 2 Marked   Age         17  57.8  5.27

datos4 %>%
    group_by(Improved) %>%
    filter(between(Age,
                   quantile(Age, 0.1),
                   quantile(Age, 0.9))) %>%
    ggboxplot(x = "Improved", y = "Age",
              add = c("mean"),
              add.params = list(color = 'red'))

Al recortar la media en un 20% vemos que se mantiene el mismo orden en las diferencias de las medias y desviación estándar para ambos grupos.

datos4 %>%
    filter(between(Age,
                   quantile(Age, 0.1),
                   quantile(Age, 0.9))) %>%
    levene_test(Age ~ Improved)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    25      1.12 0.299

Al revisar nuevamente el supuesto de homogeneidad de varianza luego de recortar la media para el uso de una prueba robusta, vemos que ahora no podemos rechazar la hipótesis nula de diferencia en la varianza y que si se cumple el supuesto.

datos4 %>%
    group_by(Improved) %>%
    filter(between(Age,
                   quantile(Age, 0.1),
                   quantile(Age, 0.9))) %>%
    shapiro_test(Age)

## # A tibble: 2 × 4
##   Improved variable statistic      p
##   <ord>    <chr>        <dbl>  <dbl>
## 1 None     Age          0.850 0.0740
## 2 Marked   Age          0.930 0.219

datos4 %>%
    group_by(Improved) %>%
    filter(between(Age,
                   quantile(Age, 0.1),
                   quantile(Age, 0.9))) %>%
    ggqqplot(x = "Age",
             facet.by = "Improved")

Luego de recortar las medias, en la prueba Shapiro-Wilk de normalidad no se rechaza la hipótesis nula de no-normalidad para ambas variables, aunque todavía en el gráfico quantil-quantil no se reflejan distribuciones normales a lo largo de la linea.

gghistogram(datos4, x = "Age",
             add = "median", rug = TRUE, bins = 15,
             color = "Improved", fill = "Improved", palette = "Dark2")

Realizando un histograma superpuesto con ambos grupos para ver la simetría de la edad, ambas distribuciones son asimétricas. Tampoco se puede utilizar la prueba U de Mann-Whitney para 2 muestras independientes.

En este escenario se considera que la prueba que estuvo más cerca de acuerdo con los supuestos fue la prueba robusta de Yuen para 2 muestras independientes, por lo que se procederá a realizar la prueba tomando en cuenta que no se rechazó la hipótesis nula de no-normalidad en la prueba Shapiro-Wilk pero que se observó que no se encontraron distribuciones normales en el gráfico quantil-quantil para ambas variables.

library(DescTools)

YuenTTest(Age ~ Improved, data = datos4)

## 
##  Yuen Two Sample t-test
## 
## data:  Age by Improved
## t = -1.1907, df = 8.4527, trim = 0.2000, p-value = 0.2662
## alternative hypothesis: true difference in trimmed means is not equal to 0
## 95 percent confidence interval:
##  -20.980775   6.604706
## sample estimates:
##   trimmed mean in group None trimmed mean in group Marked 
##                     50.88889                     58.07692

Al ver los resultados de la prueba obtenemos un p-valor mayor que 0.05, por lo que no podemos rechazar la hipótesis nula de igualdad de la media entre los grupos.

Ejercicio 5

Utiliza los datos “immer”, del paquete “MASS”, sobre el rendimiento de la cebada en los años 1931 y 1932 en un mismo campo de recolección. Comprueba mediante pruebas paramétricas, no paramétricas y robustas si han cambiado los valores medios de la producción de cebada. Interpreta y compara los resultados.

library(MASS)

data("immer")

immer %>%
     dplyr::select(Y1, Y2) %>%
     get_summary_stats(type = "mean_sd")

## # A tibble: 2 × 4
##   variable     n  mean    sd
##   <fct>    <dbl> <dbl> <dbl>
## 1 Y1          30 109.   28.7
## 2 Y2          30  93.1  24.3

library(tidyr)

immer.5 <- immer %>%
     pivot_longer(c(Y1, Y2),
                  names_to = "Year",
                  values_to = "Produccion") %>%
     arrange(Year)

ggboxplot(x = "Year", y = "Produccion",
           data = immer.5,
           add = c("mean"),
           add.params = list(color = 'red'))

immer <- immer %>% mutate(diferencias = Y1 - Y2)

immer %>% identify_outliers(diferencias)

## [1] Loc         Var         Y1          Y2          diferencias is.outlier 
## [7] is.extreme 
## <0 rows> (or 0-length row.names)

shapiro_test(immer$diferencias)

## # A tibble: 1 × 3
##   variable          statistic p.value
##   <chr>                 <dbl>   <dbl>
## 1 immer$diferencias     0.938  0.0796

ggqqplot(immer, "diferencias")

luego de ver los datos, los objetivos y analizar los supuestos, vemos que se trata de la comparación de puntuaciones para 2 muestras relacionadas, ya que el objetivo es comparar el volumen de cebada que produjeron las mismas plantaciones en dos diferentes años. Considerando que el único supuesto que no se cumplió fue el de normalidad, se determina que la prueba de hipótesis a utilizar es la prueba no paramétrica de Wilcoxon para 2 muestras relacionadas.

immer %>%
     dplyr::select(Y1, Y2) %>%
     get_summary_stats(type = "median_iqr")

## # A tibble: 2 × 4
##   variable     n median   iqr
##   <fct>    <dbl>  <dbl> <dbl>
## 1 Y1          30  103.   37.2
## 2 Y2          30   93.0  30.5

ggboxplot(x = "Year", y = "Produccion",
           data = immer.5,
           add = c("median"),
           add.params = list(color = 'red'))

immer %>% gghistogram(x = "diferencias", fill = "steelblue", add_density = TRUE)

## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.

Cuando vemos un histograma de los volúmenes registrados, la distribución tiene cierto nivel de asimetría, pero continuaremos de igual forma a ajustar el modelo.

wilcox_test(Produccion ~ Year, data = immer.5, paired = TRUE)

## # A tibble: 1 × 7
##   .y.        group1 group2    n1    n2 statistic       p
## * <chr>      <chr>  <chr>  <int> <int>     <dbl>   <dbl>
## 1 Produccion Y1     Y2        30    30      368. 0.00532

wilcox_effsize(Produccion ~ Year, data = immer.5, paired = TRUE)

## # A tibble: 1 × 7
##   .y.        group1 group2 effsize    n1    n2 magnitude
## * <chr>      <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 Produccion Y1     Y2       0.511    30    30 large

Cuando ajustamos la prueba de Wilcoxon para dos muestras relacionadas obtenemos un p-valor inferior a 0.05, por lo que podemos rechazar la hipótesis nula de igualdad y decir que existe una diferencia estadísticamente significativa entre las producciones de cebada del año 1931 y 1932.

El tamaño del efecto fue de 0.51, por lo que se trata de un tamaño del efecto grande. Podemos concluir que han cambiado los valores de la producción de cebada entre 1931 y 1932, donde en 1932 se reflejó una disminución importante.