una - ¿Por qué no puedo obtener un valor de p menor que 2.2e-16?
valor p pdf (6)
He encontrado este problema con las pruebas t y chi-cuadrado en R, pero supongo que este problema se aplica generalmente a otras pruebas. Si lo hago:
a <- 1:10
b <- 100:110
t.test(a,b)
Obtengo: t = -64.6472, df = 18.998, p-value < 2.2e-16 . Sé por los comentarios que 2.2e-16 es el valor de .Machine$double.eps , el número de punto flotante más pequeño, tal que 1 + x != 1 , pero, por supuesto, R puede representar números mucho más pequeños. También sé por la R Preguntas frecuentes que R tiene que redondear los flotantes a 53 dígitos binarios con una precisión: R Preguntas frecuentes .
Algunas preguntas: (1) ¿Tengo razón al leer que, como 53 dígitos binarios de precisión, o los valores en R < .Machine$double.eps no se calculan con precisión? (2) ¿Por qué, al realizar tales cálculos, R no proporciona un medio para mostrar un valor menor para el valor p, incluso con cierta pérdida de precisión? (3) ¿Hay una manera de mostrar un valor p más pequeño, incluso si pierdo algo de precisión? Para una sola prueba, 2 cifras decimales significativas estarían bien, para los valores que voy a corregir Bonferroni necesitaré más. Cuando digo "perder algo de precisión", creo que hay <53 dígitos binarios, pero (4) ¿estoy completamente equivocado y cualquier valor de p < .Machine$double.eps es < .Machine$double.eps inexacto? (5) ¿Es R simplemente ser honesto y otros paquetes de estadísticas no?
En mi campo, los valores p muy pequeños son la norma, algunos ejemplos: http://www.ncbi.nlm.nih.gov/pubmed/20154341 , http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.1002215 y es por eso que quiero representar valores p tan pequeños.
Gracias por su ayuda, perdón por una pregunta tan tortuosa.
Algunos paquetes R resuelven este problema. La mejor manera es a través del paquete pspearman.
source("http://www.bioconductor.org/biocLite.R")
biocLite("pspearman")
library("pspearman")
a=c(1:110,110)
b=1:111
out <- spearman.test(a, b, alternative = "greater", approximation="t-distribution")
out$p.value
[1] 3.819961e-294
Dos preguntas:
1) ¿Qué posible diferencia en la implicación estadística habría entre los valores p de 1e-16 y 1e-32? Si realmente puede justificarlo, usar los valores registrados es el camino a seguir.
2) ¿Por qué usas Wikipedia cuando estás interesado en la precisión numérica de R?
La R-FAQ dice "Otros números [que no son enteros] tienen que redondearse (por lo general) a una precisión de 53 dígitos binarios". 16 dígitos es sobre el límite. Esta es la forma de obtener los límites de precisión cuando se encuentra en la consola:
> .Machine$double.eps
[1] 2.220446e-16
Ese número es efectivamente cero cuando se interpreta en un rango de [0,1]
Estoy confundido por varias cosas en el intercambio de respuestas y comentarios aquí.
En primer lugar, cuando pruebo el ejemplo original del OP, no obtengo un valor de p tan pequeño como los que se debaten aquí (varias versiones diferentes de 2.13.x y R-devel):
a <- 1:10
b <- 10:20
t.test(a,b)
## data: a and b
## t = -6.862, df = 18.998, p-value = 1.513e-06
Segundo, cuando hago la diferencia entre los grupos mucho más grande, de hecho obtengo los resultados sugeridos por @eWizardII:
a <- 1:10
b <- 110:120
(t1 <- t.test(a,b))
# data: a and b
# t = -79.0935, df = 18.998, p-value < 2.2e-16
#
> t1$p.value
[1] 2.138461e-25
El comportamiento de la salida impresa en t.test es impulsado por su llamada a stats:::print.htest (que también es llamada por otras funciones de prueba estadística como chisq.test , como lo señala el OP), que a su vez llama format.pval , que presenta valores de p menores que su valor de eps (que es .Machine$double.eps por defecto) como < eps . Me sorprende encontrarme en desacuerdo con comentaristas tan astutos en general ...
Finalmente, aunque parece tonto preocuparse por el valor preciso de un valor p muy pequeño, el OP es correcto en el sentido de que estos valores se usan a menudo como índices de fortaleza de la evidencia en la literatura de bioinformática; por ejemplo, uno podría probar 100,000 genes candidatos y observe la distribución de los valores p resultantes (busque "gráfico de volcanes" para un ejemplo de este tipo de procedimiento).
La página de Wikipedia a la que se vinculó era para el tipo de Decimal64 que R no usa, usa dobles de edición estándar.
Primero, algunas definiciones de la página de ayuda de .Machine .
double.eps: el número de punto flotante positivo más pequeño ''x'' tal que ''1 + x! = 1''. ... Normalmente ''2.220446e-16''.
double.xmin: el número de punto flotante normalizado distinto de cero más pequeño ... Normalmente ''2.225074e-308''.
Por lo tanto, puede representar números más pequeños que 2.2e-16, pero su precisión se ve atenuada y causa problemas con los cálculos. Pruebe algunos ejemplos con números cercanos al valor representable más pequeño.
2e-350 - 1e-350
sqrt(1e-350)
Usted mencionó en un comentario que quería hacer correcciones de bonferroni. En lugar de lanzar su propio código para esto, le sugiero que use p.adjust(your_p_value, method = "bonferroni") lugar. pairwise.t.test usa esto.
Pruebe algo como este t.test(a,b)$p.value ver si eso le da la precisión que necesita. Creo que tiene más que ver con la impresión del resultado que con el valor real almacenado de la computadora que debería tener la precisión necesaria.
Recientemente tuvo el mismo problema. Fellow estadístico recomienda:
A <- cor.test(…)
p <- 2* pt(A$statistic, df = A$parameter, lower.tail=FALSE)