Valor de p NaN ao usar o goodfit de R em dados binomiais

Estou tentando testar a qualidade do ajuste de um vetor de dados de contagem para um binômio. Para fazer isso, estou usando a goodfit()função no vcdpacote. Quando executo a função, no entanto, ela retorna NaNpara o valor-p do teste Qui-quadrado. Na minha configuração, tenho um vetor de dados de contagem com 75 elementos.

> library(vcd)
> counts <- c(32, 35, 44, 35, 41, 33, 42, 49, 36, 41, 42, 45, 38, 43, 36, 
35, 40, 40, 43, 34, 39, 31, 40, 39, 36, 37, 37, 37, 32, 48, 41, 
32, 37, 36, 49, 37, 41, 36, 34, 37, 41, 32, 36, 36, 30, 33, 33, 
42, 39, 36, 36, 29, 31, 41, 36, 39, 40, 37, 39, 39, 31, 39, 37, 
40, 33, 41, 34, 46, 35, 41, 44, 38, 44, 34, 42)
> test.gof <- goodfit(counts, type="binomial", 
+                     par=list(size=length(counts), prob=0.5))

Tudo funciona bem, mas quando inspeciono o goodfit()objeto, recebo o seguinte:

> summary(test.gof)

 Goodness-of-fit test for binomial distribution

                      X^2 df  P(> X^2)
Pearson               NaN 75       NaN
Likelihood Ratio 21.48322 19 0.3107244
Warning message:
In summary.goodfit(test.gof) : Chi-squared approximation may be incorrect

Suspeitei que fosse um problema pequeno no tamanho da amostra a princípio, mas também tenho um conjunto de dados com 50 observações que não retornam NaNpara o valor p. Eu também tentei mudar o método goodfit()para ML com resultados semelhantes.

Por que essa função estaria produzindo NaNneste caso? Existe uma função alternativa para calcular o GOF em dados de contagem?

Você tem zero frequências nas contagens observadas. Isso explica os NaNseus dados. Se você olhar para o test.gofobjeto, verá que:

table(test.gof$observed)

 0  1  2  3  4  5  7  8 10 
56  5  3  2  5  1  1  2  1

você tem 56 zeros. Enfim, IMHO esta pergunta é para http://stats.stackexchange.com .

Você ficaria mais feliz com um objeto goodfit cirurgicamente alterado?

> idx <- which(test.gof$observed != 0)
> idx
 [1] 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 49 50
> test.gof$par$size <- length(  idx-1)
> test.gof$fitted <- test.gof$fitted[idx]
> test.gof$count <- test.gof$count[idx]
> test.gof$observed <- test.gof$observed[idx]
> summary(test.gof)

     Goodness-of-fit test for binomial distribution

                      X^2 df  P(> X^2)
Pearson               Inf 75 0.0000000
Likelihood Ratio 21.48322 19 0.3107244
Warning message:
In summary.goodfit(test.gof) : Chi-squared approximation may be incorrect

Tente plotá-lo. Você terá uma idéia melhor do que está acontecendo. Como mencionado anteriormente, você está recebendo NaN porque está passando 0 frequências para chisq.test ()

test.gof <- goodfit(counts, type="binomial", par=list(size=length(counts), prob=0.5)) 
plot(test.gof)
## doesn't look so good 
test.gof <- goodfit(counts, type="binomial", par=list(size=length(counts))) 
plot(test.gof)
## looks a little more clear