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
<- read_csv(url("https://zenodo.org/record/5053036/files/Winzell2004_covar.csv?download=1"))
df
<- read_csv(url("https://zenodo.org/record/5053060/files/DHS_Peru.csv?download=1")) peru
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.
<- df %>%
tab1 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()
F | M | Total | |
---|---|---|---|
(N=425) | (N=421) | (N=846) | |
Diet | |||
- chow | 225 (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 es50%
. Se puede especificar la probabilidad a ser analizada con el argumentop
.
%>%
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
<- df %>%
chi2_summ 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 paquetelsr
Independencia
Primero creamos la tabla de contingencia que vamos a analizar
<- df %>%
tab1 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()
<- tab1 %>%
chi2_summ 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