Cómo usar la prueba de ji cuadrado para determinar si los datos siguen la distribución de Poisson

La figura a continuación (Figura 1 de la página 646 de este artículo ) compara los valores observados con los valores esperados bajo la distribución de Poisson. Luego ejecuta una prueba de ji cuadrado para ver si los valores observados difieren de los valores esperados en la distribución de Poisson.

Usando R, ¿cómo es posible generar valores esperados bajo la distribución de Poisson y comparar los valores observados usando una prueba de ji cuadrado?

EDITAR:

Aquí está mi intento de hacer lo que hicieron en papel. Quiero saber si la distribución observada de variabledifiere de una distribución de Poisson. También quiero saber si lo que he hecho a continuación es el mismo procedimiento que lo que hicieron en papel. Como el valor P es> 0.05, he concluido a continuación que la distribución de variablesigue una distribución de Poisson, ¿alguien podría confirmar esto?

df <- data.frame(variable = 0:5, frequency = c(20, 10, 5, 3, 2, 1))

# estimate lambda
mean_df_variable <- mean(df$variable)

# calculate expected values if df$frequency follows a poisson distribution
library(plyr)
expected <- laply(0:5, function(x) dpois(x=x, lambda=mean_df_variable, log = FALSE))

# calculate actual distribution of df$frequency
observed <- df$frequency/sum(df$frequency)

# does distribution of df$frequency differ from a poisson distribution? Apparently 
#   not because P-value is > 0.05
chisq.test(expected, observed)

La forma en que realizó la prueba de chi-cuadrado no es correcta. Hay varios problemas. Primero, su marco de datos se ve así:

  variable frequency
1        0        20
2        1        10
3        2         5
4        3         3
5        4         2
6        5         1

Entonces, cuando corres mean(df$variable), obtienes 2.5, que es solo la media de 0:5. Es decir, no tiene peso. En cambio, cree su variable así:

x = rep(0:5, times=c(20, 10, 5, 3, 2, 1))
table(x)
# x
#  0  1  2  3  4  5 
# 20 10  5  3  2  1
mean(x)
# [1] 1.02439

La table()llamada muestra que el código nos da lo que queríamos, por lo que mean()estima lambda correctamente.

A continuación, sus probabilidades estimadas solo van a 5, pero la distribución de Poisson va al infinito. Por lo tanto, debe tener en cuenta las probabilidades de los valores que no tiene en su conjunto de datos. Esto no es difícil de hacer, solo calcula el complemento:

probs = dpois(0:5, lambda=mean(x))
probs
# [1] 0.359015310 0.367771781 0.188370912 0.064321775 0.016472650 0.003374884
comp = 1-sum(probs)
# [1] 0.0006726867

Por último, en Rla chisq.test()función de, los argumentos x=y y=no son exactamente para los valores esperados y observados en la forma en que configura esto. Por un lado, lo que llama "esperado" son en realidad probabilidades (es decir, la salida de dpois()), para hacer estos valores esperados, tendría que multiplicar esas probabilidades (y asegurarse de incluir el cumplido) por el recuento total. Pero incluso entonces, no los usarías para y=. En cualquier caso, en realidad no tiene que hacer eso, solo puede asignar las probabilidades al p=argumento. Además, deberá agregar un 0vector de valores observados para representar todos los valores posibles que no aparecen en su conjunto de datos:

chisq.test(x=c(20, 10, 5, 3, 2, 1, 0), p=c(probs, comp))

#  Chi-squared test for given probabilities
# 
# data:  c(20, 10, 5, 3, 2, 1, 0)
# X-squared = 12.6058, df = 6, p-value = 0.04974
# 
# Warning message:
#   In chisq.test(x = c(20, 10, 5, 3, 2, 1, 0), p = c(probs, comp)) :
#   Chi-squared approximation may be incorrect

El mensaje de advertencia sugiere que preferimos simular en su lugar, así que lo intentamos nuevamente:

chisq.test(x=c(20, 10, 5, 3, 2, 1, 0), p=c(probs, comp), simulate.p.value=TRUE)

# Chi-squared test for given probabilities with simulated p-value 
#   (based on 2000 replicates)
# 
# data:  c(20, 10, 5, 3, 2, 1, 0)
# X-squared = 12.6058, df = NA, p-value = 0.07046

Presumiblemente, este es un valor p más preciso, pero plantea una pregunta sobre cómo debe interpretarse. Usted pregunta "Como el valor P es> 0.05, he concluido a continuación que la distribución de la variable sigue una distribución de Poisson, ¿alguien podría confirmar esto?" Usando el enfoque correcto, notamos que el primer valor p fue solo <.05, pero el segundo valor p (simulado) fue solo> .05. Aunque el último valor p es más preciso, no me apresuraría a concluir que los datos provienen de una distribución de Poisson. Aquí hay algunos hechos a tener en cuenta:

• Como se indica en el título de un artículo de Gelman y Stern, la diferencia entre "significativo" y "no significativo" no es estadísticamente significativa .
• Los datos reales nunca provienen de ninguna de las distribuciones idealizadas, un hecho señalado por @Glen_b recientemente aquí: ¿A qué distribución se parece este histograma?

  Los datos reales tienden a no seguir las formas distributivas simples de las distribuciones comunes de uno, dos o tres parámetros. Las distribuciones reales son más como mezclas heterogéneas. Las formas distributivas simples son ficciones convenientes (modelos, para ser precisos): se aproximan a la realidad de una manera que facilita su trabajo.

• No puede utilizar el hecho de un resultado no significativo para afirmar la hipótesis nula, como explico aquí: ¿Por qué los estadísticos dicen que un resultado no significativo significa "no se puede rechazar el nulo" en lugar de aceptar la hipótesis nula?

Si he entendido lo que querías decir, deberías:

1. estimar el parámetro de la distribución de Poisson para sus datos, suponiendo que esté distribuido por Poisson, digamos

lambdaEst = mean(x)

2. calcular, por cada , sus probabilidades teóricas suponiendo una distribución de Poisson, por ejemplo$0,1, 2, ...$

probTheo0 = dpois(x = 0, lambda = lambdaEst, log = FALSE)

3. luego compare las probabilidades reales con las teóricas a través de una prueba de chi-cuadrado siguiendo este enfoque ChiSquare Test CV solution