在适合卡方检验的优度上获得 p 值 = 1

Question

我正在尝试使用 R 对一系列观察结果进行泊松拟合优度检验。我正在计算每分钟有多少人在 57 分钟内做了某件事。我从来没有得到任何大于 13 的观察值，我得到了以下数据： （适用于 0 至 13 人以上的情况）：

observed = c(3/57, 4/57, 9/57, 7/57, 9/57, 8/57, 2/57, 3/57, 7/57, 2/57, 1/57, 0, 1/57, 1/57, 0)

表示我观察了 3 次 0 人、4 次 1 人、9 次 2 人等等（最后的 0 表示我从未见过 14 人或更多人）。

mn = 4.578947 
cases = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)
estimated = c()
for (i in cases)(estimated <- c(estimated, dpois(i, lambda = mn)))
estimated <- c(estimated, (1-ppois(13, lambda=mn)))

其中mn 是从数据中获得的平均值。 最后，我跑了

 chisq.test(observed, p=estimated)

我得到：

 Chi-squared test for given probabilities

data:  observed
X-squared = 1.0182, df = 14, p-value = 1

Warning message:
In chisq.test(observed, p = estimated) :
  Chi-squared approximation may be incorrect

我在这方面并不精通（既不是统计数据，也不是 R 编程），但我认为我不应该得到正好为 1.0 的 p 值。我究竟做错了什么？ （顺便说一句：我的代码很可能不是我想要做的最佳选择，但我几乎不使用 R 并且它不是我现在工作的重点。）

Answer 1

您的观察值应该是计数，而不是比例：

> chisq.test(observed*57, p=estimated)

    Chi-squared test for given probabilities

data:  observed * 57
X-squared = 58.036, df = 14, p-value = 2.585e-07

根据chisq.test 的 R 帮助文件：

如果 x 是一个包含一行或一列的矩阵，或者如果 x 是一个向量并且 y 是 未给出，则执行拟合优度检验（x 被视为 一维列联表）。 x 的条目必须是 非负整数。

（强调我的）

您可以使用手册中的一些示例代码对此进行测试

应该怎么做：

> x <- c(89,37,30,28,2)
> p <- c(0.40,0.20,0.20,0.19,0.01)
> chisq.test(x, p = p)

    Chi-squared test for given probabilities

data:  x
X-squared = 5.7947, df = 4, p-value = 0.215

Warning message:
In chisq.test(x, p = p) : Chi-squared approximation may be incorrect

并且犯了和你一样的错误：

> chisq.test(x/sum(x), p = p)

    Chi-squared test for given probabilities

data:  x/186
X-squared = 0.031154, df = 4, p-value = 0.9999

Warning message:
In chisq.test(x/186, p = p) : Chi-squared approximation may be incorrect

Answer 2

首先，进行拟合优度测试，观察频率和 bin概率是必需的。

 observed = c(3, 4, 9, 7, 9, 8, 2, 3, 7, 2, 1, 0, 1, 1, 0)       # keep counts

概率是正确的：

 mn = 4.578947 
 prob = c()
 for (i in cases)     (prob <- c(prob, dpois(i, lambda = mn)))
 prob <- c(prob, (1-ppois(13, lambda=mn)))           # prob for 13 and plus category

最重要的是，分类/类别中的预期频率至少应为 5。 Chisq-test 对小样本无效。 这就是为什么您会收到警告 （请参阅类别 1,2 和 8-15 的预期频率）：

poisson_df <- data.frame(observed, prob)
poisson_df$expected = sum(poisson_df$observed)*poisson_df$prob

poisson_df

#   observed   prob          expected
#1         3   0.0102657004  0.58514492
#2         4   0.0470060980  2.67934759
#3         9   0.1076192157  6.13429530
#4         7   0.1642608950  9.36287101
#5         9   0.1880354831 10.71802253
#6         8   0.1722009022  9.81545143
#7         2   0.1314164674  7.49073864
#8         3   0.0859641485  4.89995646
#9         7   0.0492031600  2.80458012
#10        2   0.0250331846  1.42689152
#11        1   0.0114625626  0.65336607
#12        0   0.0047714970  0.27197533
#13        1   0.0018207026  0.10378005
#14        1   0.0006413001  0.03655410
#15        0   0.0002986829  0.01702492

chisq.test(x = poisson_df$observed, p= poisson_df$prob)

# Chi-squared test for given probabilities

# data:  observed
# X-squared = 58.036, df = 14, p-value = 2.585e-07

Warning message:
In chisq.test(x = poisson_df$observed, p= poisson_df$prob) :
Chi-squared approximation may be incorrect

因此，需要适当地创建分箱。需要注意的是Chisq-test对分箱很敏感em>，bin的一种方式如下：

cat_eq_3_less <- apply(poisson_df[1:3,], 2 , sum)        # sum of 1 to 3 categories
cat_eq_8_plus <- apply(poisson_df[8:15,], 2 , sum)       # sum 8 to 15 categories

corrected_df <- rbind(cat_eq_3_less, poisson_df[4:7,], cat_eq_8_plus)

 corrected_df
 #   observed     prob       expected
 #        16      0.1648910  9.398788
 #         7      0.1642609  9.362871
 #         9      0.1880355 10.718023
 #         8      0.1722009  9.815451
 #         2      0.1314165  7.490739
 #        15      0.1791952 10.214129

chisq.test(x = corrected_df$observed, p = corrected_df$prob)

Chi-squared test for given probabilities

data:  corrected_df$observed
X-squared = 12.111, df = 5, p-value = 0.0333