Prueba de hipótesis para datos categóricos: Chi cuadrado. Test exacto de Fisher - FELSOCEM

Gabriel Carrasco-Escobar

Paquetes y Data

Los paquetes que se utilizaran son:

library(tidyverse)
library(janitor)
library(pubh)
library(ggmosaic)
library(rstatix)
library(Publish)

Para los ejercicios a continuacion se utilizaron las siguientes bases de datos:

  • Winzell2004_covar.csv
  • DHS_Peru.csv
df <- read_csv(url("https://zenodo.org/record/5053036/files/Winzell2004_covar.csv?download=1"))

peru <- read_csv(url("https://zenodo.org/record/5053060/files/DHS_Peru.csv?download=1"))

Medidas de Resumen

Tablas de frecuencias

El paquete janitor tiene la función tabyl que permite generar tablas de frecuencias y se integra en el ecosistema del tidyverse.

df %>%
  tabyl(Diet)
##  Diet   n   percent
##  chow 449 0.5307329
##    hf 397 0.4692671

tabyl permite almacenar el objeto del cálculo de frecuencias en formato tidy para su manipulación posterior.

tab1 <- df %>%
  tabyl(Diet)
Diet n percent
chow 449 0.5307329
hf 397 0.4692671

Adicionalmente podemos configurar nuestra tabla con adorn_*(). Por ejemplo, podemos utilizar adorn_totals("row") para agregar el total de las filas y adorn_pct_formatting() para dar formato a la columna de porcentajes.

df %>%
  tabyl(Diet) %>%
  adorn_totals("row") %>%
  adorn_pct_formatting()
Diet n percent
chow 449 53.1%
hf 397 46.9%
Total 846 100.0%

Tablas de contingengias

Podemos realizar tablas de doble entrada (frecuencias de 2 variables) con tabyl

df %>%
  tabyl(Diet, Sex)
Diet F M
chow 225 224
hf 200 197

Al igual que con las freuencias de una vairable podemos cambiar el formato de la tabla con adorn_*(). Más información de las opciones de formato se encuentra detallado en la documentación de janitor.

df %>%
  tabyl(Diet, Sex) %>%
  adorn_percentages("col") %>%
  adorn_pct_formatting(digits = 2) %>%
  adorn_ns()
Diet F M
chow 52.94% (225) 53.21% (224)
hf 47.06% (200) 46.79% (197)

Otro paquete con funciones útiles para medidas descriptivas es pubh (enfocado en salud pública). La función cross_tab() permite realizar también tablas de doble entrada.

df %>%
  cross_tab(Sex ~ Diet) %>%
  theme_pubh()
FMTotal
(N=425)(N=421)(N=846)
Diet
- chow225 (52.9%)224 (53.2%)449 (53.1%)
- hf 200 (47.1%)197 (46.8%)397 (46.9%)

Visualizacion y Analisis de datos Exploratorio (ADE)

Frecuencias

Usaremos la geometría geom_bar() del paquete ggplot para generar un gráfico de barras de una variable (tabulación de un solo sentido).

En el siguiente ejemplo, exploraremos la variable Diet para visualizar la distribución de las freqcuencias de dietas chow y hf.

df %>%
  ggplot(aes(x=Diet, fill = Diet)) +
  geom_bar()

Mosaicos

Los gráficos de mosaicos nos permite comparar dos variables categoricas. Usaremos la geometría geom_mosaic() del paquete ggmosaic

df %>%
  ggplot() + 
  geom_mosaic(aes(x = product(Sex, Diet), fill = Sex))

Pruebas de hipotesis

Chi2

Bondad de Ajuste (Goodness of fit)

Usaremos la prueba de Chi-cuadrado con el paquete rstatix y las funciones chisq_test() para comparar los valores observados contra una probabilidad de referencia.

Si no se especifica la probabilidad de referencia, chisq_test() utiliza un probabilidad homogenea entre las catagorias de la variable categórica. Para una variable binaria la probablidad usada es 50%. Se puede especificar la probabilidad a ser analizada con el argumento p.

df %>%
  tabyl(Diet) %>%
  pull(n) %>%
  chisq_test()
n statistic p df method p.signif
2 3.196217 0.0738 1 Chi-square test ns

Podemos explorar a profundidad la prueba de Chi-cuadrado con la función chisq_descriptives().

df %>%
  tabyl(Diet) %>%
  pull(n) %>%
  chisq_test() %>%
  chisq_descriptives() %>%
  adorn_totals("row")
Var1 observed prop expected resid std.resid
A 449 0.5307329 423 1.264163 1.787797
B 397 0.4692671 423 -1.264163 -1.787797
Total 846 1.0000000 846 0.000000 0.000000

El valor del estadístico de Pearson's Chi-squared es la suma de los residuos elevados al cuadrado

chi2_summ <- df %>%
  tabyl(Diet) %>%
  pull(n) %>%
  chisq_test() %>%
  chisq_descriptives()

sum(chi2_summ$resid^2)
## [1] 3.196217

Tambien se puede usar la función goodnessOfFitTest() del paquete lsr

Independencia

Primero creamos la tabla de contingencia que vamos a analizar

tab1 <- df %>%
  tabyl(Diet, Sex)
Diet F M
chow 225 224
hf 200 197

El paquete janitor tiene una extensión de la función chisq.test() para usar en objetos de clase tabyl.

chisq.test(tab1, cor = F)
## 
##  Pearson's Chi-squared test
## 
## data:  tab1
## X-squared = 0.0059848, df = 1, p-value = 0.9383

Podemos también usar la función chisq_test() del paquete rstatix.

tab1 %>%
  column_to_rownames("Diet") %>%
  chisq_test(cor = F)
n statistic p df method p.signif
846 0.0059848 0.938 1 Chi-square test ns

Podemos explorar a profundidad la prueba de Chi-cuadrado con el paquete rstatix y la función chisq_descriptives()

chi2_summ <- tab1 %>%
  column_to_rownames("Diet") %>%
  chisq_test(cor = F) %>%
  chisq_descriptives()
Var1 Var2 observed prop row.prop col.prop expected resid std.resid
chow F 225 0.2659574 0.5011136 0.5294118 225.5615 -0.0373844 -0.0773615
hf F 200 0.2364066 0.5037783 0.4705882 199.4385 0.0397575 0.0773615
chow M 224 0.2647754 0.4988864 0.5320665 223.4385 0.0375616 0.0773615
hf M 197 0.2328605 0.4962217 0.4679335 197.5615 -0.0399459 -0.0773615

El valor del estadístico de Pearson's Chi-squared es la suma de los residuos elevados al cuadrado

sum(chi2_summ$resid^2)
## [1] 0.005984801

Fisher

El paquete janitor tiene una extensión de la función fisher.test() para usar en objetos de clase tabyl.

fisher.test(tab1)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab1
## p-value = 0.9452
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7483337 1.3081149
## sample estimates:
## odds ratio 
##  0.9894081

Podemos también usar la función fisher_test() del paquete rstatix.

tab1 %>%
  column_to_rownames("Diet") %>%
  fisher_test()
n p p.signif
846 0.945 ns

Tablas epidemiologicas

También podemos usar la función contingency() del paquete pubh para generar una salida de datos más completa

df %>%
  contingency(Sex ~ Diet)
##          Outcome
## Predictor   M   F
##      hf   197 200
##      chow 224 225
## 
##              Outcome +    Outcome -      Total        Inc risk *        Odds
## Exposed +          197          200        397              49.6       0.985
## Exposed -          224          225        449              49.9       0.996
## Total              421          425        846              49.8       0.991
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                               0.99 (0.87, 1.14)
## Odds ratio                                   0.99 (0.76, 1.30)
## Attrib risk *                                -0.27 (-7.02, 6.48)
## Attrib risk in population *                  -0.13 (-5.85, 5.60)
## Attrib fraction in exposed (%)               -0.54 (-15.15, 12.22)
## Attrib fraction in population (%)            -0.25 (-6.82, 5.92)
## -------------------------------------------------------------------
##  Test that OR = 1: chi2(1) = 0.006 Pr>chi2 = 0.94
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units 
## 
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  dat
## X-squared = 7.1725e-05, df = 1, p-value = 0.9932

Otro paquete útil para mostrar múltiples pruebas de hipótesis es table2x2() del paquete Publish

tab1 %>%
  column_to_rownames("Diet") %>%
  table2x2()
## _____________________________
## 
## 2x2 contingency table
## _____________________________
## 
##            F       M     Sum
## chow     225     224     449
## hf       200     197     397
## --        --      --      --
## Sum      425     421     846
## 
## _____________________________
## 
## Statistics
## _____________________________
## 
## 
## a= 225 
## b= 224 
## c= 200 
## d= 197
## 
## p1=a/(a+b)= 0.5011 
## p2=c/(c+d)= 0.5038 
## 
## _____________________________
## 
## Risk difference
## _____________________________
## 
## Risk difference = RD = p1-p2 = -0.002665
## Standard error = SE.RD = sqrt(p1*(1-p1)/(a+b)+p2*(1-p2)/(c+d)) = 0.03445
## Lower 95%-confidence limit: = RD - 1.96 * SE.RD = -0.07018
## Upper 95%-confidence limit: = RD + 1.96 * SE.RD = 0.06485
## 
## The estimated risk difference is -0.3%  (CI_95%: [-7.0;6.5]).
## 
## _____________________________
## 
## Risk ratio
## _____________________________
## 
## Risk ratio = RR = p1/p2 = 0.9947
## Standard error = SE.RR = sqrt((1-p1)/a+(1-p2)/c)= 0.9947
## Lower 95%-confidence limit: = RR * exp(- 1.96 * SE.RR) = 0.8697
## Upper 95%-confidence limit: = RR * exp(1.96 * SE.RR) = 1.1377
## 
## The estimated risk ratio is 0.995 (CI_95%: [0.870;1.138]).
## 
## _____________________________
## 
## Odds ratio
## _____________________________
## 
## Odds ratio = OR = (p1/(1-p1))/(p2/(1-p2)) = 0.9894
## Standard error = SE.OR = sqrt((1/a+1/b+1/c+1/d)) = 0.1378
## Lower 95%-confidence limit: = OR * exp(- 1.96 * SE.OR) = 0.7552
## Upper 95%-confidence limit: = OR * exp(1.96 * SE.OR) = 1.2961
## 
## The estimated odds ratio is 0.989 (CI_95%: [0.755;1.296]).
## 
## _____________________________
## 
## Chi-square test
## _____________________________
## 
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table2x2
## X-squared = 7.1725e-05, df = 1, p-value = 0.9932
## 
## 
## _____________________________
## 
## Fisher's exact test
## _____________________________
## 
## 
##  Fisher's Exact Test for Count Data
## 
## data:  table2x2
## p-value = 0.9452
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7483337 1.3081149
## sample estimates:
## odds ratio 
##  0.9894081