Comment utiliser le test du chi carré pour déterminer si les données suivent la distribution de Poisson

La figure ci-dessous (figure 1 de la page 646 de cet article ) compare les valeurs observées aux valeurs attendues sous la distribution de Poisson. Il exécute ensuite un test du chi carré pour voir si les valeurs observées diffèrent des valeurs attendues sous la distribution de Poisson.

En utilisant R, comment est-il possible de générer des valeurs attendues sous distribution de Poisson et de comparer les valeurs observées à l'aide d'un test du chi carré?

ÉDITER:

Voici ma tentative de faire ce qu'ils ont fait sur papier. Je veux savoir si la distribution observée de variablediffère d'une distribution de Poisson. Je veux également savoir si ce que j'ai fait ci-dessous est la même procédure que ce qu'ils ont fait sur papier. Comme la valeur P est> 0,05, j'ai conclu ci-dessous que la distribution de variablesuit une distribution de Poisson - quelqu'un pourrait-il le confirmer?

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 façon dont vous avez fait le test du chi carré n'est pas correcte. Il y a plusieurs problèmes. Tout d'abord, votre bloc de données ressemble à ceci:

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

Donc, quand vous courez mean(df$variable), vous obtenez 2.5, ce qui est juste la moyenne de 0:5. Autrement dit, il n'est pas pondéré. Au lieu de cela, créez votre variable comme ceci:

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

L' table()appel montre que le code nous donne ce que nous voulions, et mean()estime donc lambda correctement.

Ensuite, vos probabilités estimées ne vont qu'à 5, mais la distribution de Poisson va à l'infini. Vous devez donc tenir compte des probabilités des valeurs que vous n'avez pas dans votre ensemble de données. Ce n'est pas difficile à faire, il suffit de calculer le complément:

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

Enfin, dans Rla chisq.test()fonction de, les arguments x=et y=ne correspondent pas exactement aux valeurs attendues et observées dans la façon dont vous configurez cela. D'une part, ce que vous appelez "attendu" sont en fait des probabilités (c'est-à-dire la sortie de dpois()), pour faire ces valeurs attendues, vous devrez multiplier ces probabilités (et assurez-vous d'inclure le compliment) par le nombre total. Mais même alors, vous ne les utiliseriez pas y=. Quoi qu'il en soit, vous n'avez pas vraiment à le faire, vous pouvez simplement affecter les probabilités à l' p=argument. De plus, vous devrez ajouter un 0à votre vecteur de valeurs observées pour représenter toutes les valeurs possibles qui n'apparaissent pas dans votre jeu de données:

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

Le message d'avertissement suggère que nous préférons peut-être simuler à la place, nous essayons donc à nouveau:

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

Il s'agit vraisemblablement d'une valeur de p plus précise, mais cela soulève une question sur la façon de l'interpréter. Vous demandez "Comme la valeur P est> 0,05, j'ai conclu ci-dessous que la distribution de la variable suit une distribution de Poisson - quelqu'un pourrait-il confirmer cela?" En utilisant la bonne approche, nous notons que la première valeur de p était juste <0,05, mais la deuxième valeur de p (simulée) était juste> 0,05. Bien que cette dernière valeur de p soit plus précise, je ne me précipiterais pas pour conclure que les données provenaient d'une distribution de Poisson. Voici quelques faits à garder à l'esprit:

• Comme indiqué dans le titre d'un article de Gelman et Stern, la différence entre "significatif" et "non significatif" n'est pas elle-même statistiquement significative .
• Les données réelles ne proviennent jamais d'aucune des distributions idéalisées, un fait bien souligné par @Glen_b récemment ici: à quelle distribution cet histogramme ressemble-t-il?

  Les données réelles ont tendance à ne pas suivre les formes distributionnelles simples des distributions communes à un, deux ou trois paramètres. Les distributions réelles ressemblent davantage à des mélanges hétérogènes. Les formes de distribution simples sont des fictions pratiques (modèles, pour être précis) - elles se rapprochent de la réalité de manière à en faciliter l'utilisation.

• Vous ne pouvez pas utiliser le fait d'un résultat non significatif pour affirmer l'hypothèse nulle, comme je l'explique ici: pourquoi les statisticiens disent-ils qu'un résultat non significatif signifie "vous ne pouvez pas rejeter le nul" au lieu d'accepter l'hypothèse nulle?

Si j'ai compris ce que vous vouliez dire, vous devriez:

1. estimer le paramètre de la distribution de Poisson pour vos données, en supposant qu'il est distribué par Poisson, disons

lambdaEst = mean(x)

2. calculer, pour chaque , leurs probabilités théoriques en supposant une distribution de Poisson, par exemple$0,1, 2, ...$

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

3. puis comparer les probabilités réelles et théoriques via un test du chi carré suivant cette approche ChiSquare Test CV solution