¿Cómo calcular el valor p de una razón de probabilidades en R?

Tengo la siguiente tabla de valores:

25  75
38  162

La razón de posibilidades es 0.7037 y el registro (OR) es -0.3514. Para una tabla de contingencia con valores a, b, c y d, la varianza de log (OR) viene dada por

(1/a + 1/b + 1/c + 1/d)

¿Cómo puedo calcular el p.value de log (OR) a partir de estos datos en R (si es significativamente diferente de 0)?

r variance rnso
fuente

Respuestas:

Puede usar la prueba exacta de Fisher, que ingresa una tabla de contingencia y genera un valor p, con una hipótesis nula de que la razón de probabilidades es 1 y una hipótesis alternativa de que la razón de probabilidades no es igual a 1.

(tab <- matrix(c(38, 25, 162, 75), nrow=2))
#      [,1] [,2]
# [1,]   38  162
# [2,]   25   75
fisher.test(tab)
# 
#   Fisher's Exact Test for Count Data
# 
# data:  tab
# p-value = 0.2329
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  0.3827433 1.3116294
# sample estimates:
# odds ratio 
#  0.7045301

En este caso el valor p es 0.23.

josliber
fuente

Gracias por una forma inteligente de determinar el valor p. La prueba de chi-cuadrado también se puede usar de manera similar.

rnso

@rnso seguro, aunque la prueba exacta de Fisher se prefiere sobre Chi-cuadrado cuando tiene tamaños de celda pequeños en su tabla de contingencia.

josliber

Este es un mito antiguo, pero desafortunadamente no es cierto. El ordinario Pearson

χ^{2}

$\chi^2$ proporciona más precisa

P

$P$ -valores que la llamada prueba "exacta" de Fisher incluso cuando las frecuencias esperadas son tan bajas como 1.0.

Frank Harrell

¿podrías decir un poco más sobre esto @FrankHarrell? Sé que el

χ^{2}

$\chi^2$ sería un resultado asintótico, mientras que la prueba exacta de Fisher se basa en la distribución exacta, ¿cómo es el

p

$p$ -valor más "preciso" utilizando el método asintótico?

bdeonovic

Vea comentarios extensos sobre esto en el sitio. Brevemente, los valores P de la prueba de Fisher son demasiado grandes. El error absoluto medio en los valores P de la prueba de Pearson es menor. Fisher es solo "exacto" en el sentido de que los valores P están "garantizados" para no ser demasiado pequeños.

Frank Harrell

Otra forma de hacerlo (que no sea la prueba exacta de Fisher) es poner los valores en un GLM binomial:

d <- data.frame(g=factor(1:2),
                s=c(25,75),
                f=c(38,162))
g <- glm(s/(s+f)~g,weights=s+f,data=d,
    family="binomial")
coef(summary(g))["g2",c("Estimate","Pr(>|z|)")]
##   Estimate   Pr(>|z|) 
## -0.3513979  0.2303337

Para obtener la prueba de razón de probabilidad (un poco más precisa que la de Wald $p$ -valor mostrado arriba), hacer

anova(g,test="Chisq")

lo que da

##      Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                     1     1.4178         
## g     1   1.4178         0     0.0000   0.2338

(LRT $p=0.2338 \approx$ Wald $p=0.2303337 \approx$ Pescador $p=0.2329$ en este caso porque la muestra es bastante grande)

Ben Bolker
fuente

Es mejor generalizar la solución y usar la razón de probabilidad $\chi^2$ prueba de un modelo estadístico como el modelo logístico. La prueba LR proporciona bastante precisa $P$ -valores. Esto también maneja casos en los que necesita probar más de un parámetro, por ejemplo, problemas de 3 grupos, efectos continuos que no son lineales, etc. La prueba LR para el modelo general (que es todo lo que se necesita en este ejemplo ya que no hay ajuste variables) pueden obtenerse fácilmente en la base R o utilizando el rmspaquete, p. ej.

f <- lrm(y ~ groups, weights=freqs)
f  # prints LR chi-sq, d.f., P, many other quantities

Aquí los modelos anidados son este modelo y un modelo de solo intercepción.

Frank Harrell
fuente

Podría encontrar que la prueba LR (lrtest) se usa para comparar modelos anidados. ¿Cómo podemos usarlo aquí? ¿Podría escribir una línea de código R para ello?

rnso

por lo que vale, este es más o menos el mismo enfoque estadístico (aunque con una mejor explicación) que en mi respuesta anterior. lrm()tiene diferentes valores predeterminados, formatos de salida, etc., pero el modelo estadístico (IIUC) es el mismo queglm(...,family="binomial")

Ben Bolker