Prueba exacta de Fisher o regresión logística para determinar si las distribuciones difieren

Question

Tengo los siguientes datos:

                       Somkers
                     Yes   No  All
Groups   ]10,20]      35   6   41
         ]20,30]      20  13   33
         ]30,40]      10  15   25
         ]40,50]      15   9   24
         All          80  43   123

Entiendo que este tipo de datos es una tabla de 4 por 2 y aquí se puede aplicar una regresión logística. Mi objetivo aquí es probar si los fumadores y los no fumadores tienen la misma distribución en los 4 grupos. ¿Qué prueba es más adecuada en este caso: una prueba de homogeneidad o la prueba exacta de Fisher? ¿Tengo que especificar completamente el modelo logístico para poder realizar estas pruebas?

Answer 1

La prueba chi-cuadrado será la más sencilla y adecuada. La prueba exacta de Fisher evalúa las diferencias condicionadas a márgenes fijos, lo que casi con toda seguridad es inadecuado en este caso. La regresión logística estaría bien, pero la chi-cuadrado sería más sencilla; además, la LR evalúa realmente el tabaquismo en función de los grupos, lo que no se ajusta conceptualmente a la pregunta real.

d = read.table(text="Groups     Yes   No  All
                    ']10,20]'    35    6   41
                    ']20,30]'    20   13   33
                    ']30,40]'    10   15   25
                    ']40,50]'    15    9   24", header=T)

tab = as.table(as.matrix(d[,-c(1,4)]))
names(dimnames(tab)) = c("Groups", "Smoker")
rownames(tab)        = d[,1]
colnames(tab)        = names(d)[2:3]

chisq.test(tab)
#   Pearson's Chi-squared test
# 
# data:  tab
# X-squared = 14.697, df = 3, p-value = 0.002095

Permítanme añadir un par de notas adicionales: La prueba chi-cuadrado sólo le da un valor p para el nulo de que las distribuciones son las mismas. Es posible que desee caracterizar cómo difieren. Un par de maneras de hacer esto sería hacer una tabla de proporciones por columnas, y hacer un gráfico de mosaico:

round(prop.table(tab, 2), 3)
#          Smoker
# Groups      Yes    No
#   ]10,20] 0.438 0.140
#   ]20,30] 0.250 0.302
#   ]30,40] 0.125 0.349
#   ]40,50] 0.188 0.209
windows()
  mosaicplot(t(tab), shade=T)

De ellas se desprende que hay "muy pocos" no fumadores en la ]10,20] y "demasiados" en el grupo ]30,40] grupo.

---

Si quisiera utilizar una regresión logística para probar estos datos. Es muy sencillo:

mod = glm(cbind(Yes, No)~Groups, d, family=binomial)
summary(mod)
# Call:
# glm(formula = cbind(Yes, No) ~ Groups, family = binomial, data = d)
# 
# Deviance Residuals: 
# [1]  0  0  0  0
# 
# Coefficients:
#               Estimate Std. Error z value Pr(>|z|)    
# (Intercept)     1.7636     0.4419   3.991 6.57e-05 ***
# Groups]20,30]  -1.3328     0.5676  -2.348 0.018866 *  
# Groups]30,40]  -2.1691     0.6016  -3.606 0.000311 ***
# Groups]40,50]  -1.2528     0.6108  -2.051 0.040249 *  
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 1.5415e+01  on 3  degrees of freedom
# Residual deviance: 2.4425e-15  on 0  degrees of freedom
# AIC: 22.657
# 
# Number of Fisher Scoring iterations: 3

La prueba de Groups no coincide con ninguna de las pruebas individuales mostradas en la salida de resumen. Con una sola variable en el modelo, la prueba de la variable en su conjunto es la misma que la prueba del modelo en su conjunto. Desafortunadamente, R no le ofrece esto por defecto como lo hace para un modelo lineal. Sin embargo, puede utilizar las desviaciones nula y residual y los grados de libertad para obtener una prueba de razón de verosimilitud:

1-pchisq(mod$null.deviance-deviance(mod), df=mod$df.null-mod$df.residual)
# [1] 0.001494036

Con una sola variable en el modelo, una forma conveniente de obtener la prueba del modelo en su conjunto es utilizar anova.glm() . Estableciendo test="LRT" se obtiene lo mismo que con el método manual anterior, y utilizando test="Rao" , se obtiene la prueba de puntuación, que es la misma que la chi-cuadrado en la parte superior:

anova(mod, test="LRT")
# Analysis of Deviance Table
# Model: binomial, link: logit
# Response: cbind(Yes, No)
# Terms added sequentially (first to last)
# 
#        Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
# NULL                       3     15.415            
# Groups  3   15.415         0      0.000 0.001494 **
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(mod, test="Rao")
# Analysis of Deviance Table
# Model: binomial, link: logit
# Response: cbind(Yes, No)
# Terms added sequentially (first to last)
# 
#        Df Deviance Resid. Df Resid. Dev    Rao Pr(>Chi)   
# NULL                       3     15.415                   
# Groups  3   15.415         0      0.000 14.697 0.002095 **
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

Por último, la denominación de sus filas (es decir, el Groups ) es sospechoso. ¿Es el resultado de categorizar una variable originalmente continua? Si es así, es muy no recomendada. La categorización puede crear arbitrariamente la apariencia de diferentes distribuciones cuando en realidad no existen. Para hacerse una idea, puede ser útil leer esta excelente respuesta (aunque el contexto difiere): Evaluar la distribución aproximada de los datos a partir de un histograma . Sería mucho mejor utilizar los valores continuos originales y decidir a priori qué tipo de diferencia en las distribuciones le interesaría detectar (desplazamiento de la media, diferencia en la dispersión, asimetría, colas gruesas, comportamiento de otras colas, multimodalidad, etc.), y comprobarlo explícitamente.