11
$\begingroup$

I am attempting to test the goodness of fit for a vector of count data to a binomial. To do so I am using the goodfit() function in the vcd package. When I run the function, however, it returns NaN for the p-value of the Chi-squared test. In my setup, I have a vector of count data with 75 elements.

> 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))

Everything works fine, but when I inspect the goodfit() object I get the following:

> 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

I suspected it was a small sample size issue at first, but I also have a data set with 50 observations that does not return NaN for the p-value. I have also tried to switch the method in goodfit() to ML with similar results.

Why would this function be producing NaN in this case? Is there an alternative function to calculate GOF on count data?

$\endgroup$
4
  • $\begingroup$ @Gavin: should it be closed or could it be moved? $\endgroup$ Apr 27, 2011 at 18:44
  • $\begingroup$ @Joshua If it can be moved by a mod, moved would be best as @DrewConway has done a reasonable job with the Q and reproducibility. How do we achieve a migrate? Flag the post? $\endgroup$ Apr 27, 2011 at 18:48
  • $\begingroup$ Damn, I was writing an answer and had to move for 10 mins from machine. And now I came back and saw your comments. =) $\endgroup$
    – aL3xa
    Apr 27, 2011 at 18:56
  • $\begingroup$ This question has been merged with its near-duplicate. $\endgroup$
    – whuber
    Apr 28, 2011 at 2:10

3 Answers 3

7
$\begingroup$

You have zero frequencies in observed counts. That explains NaNs in your data. If you look at test.gof object, you'll see that:

table(test.gof$observed)

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

you have 56 zeros. Anyway, IMHO this question is for http://stats.stackexchange.com.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks for the re-post from SO, but how do you recommend I deal with the zeroes given that I still want to estimate a GOF? $\endgroup$
    – DrewConway
    Apr 28, 2011 at 13:04
  • $\begingroup$ Hi Drew, sorry for my late answer. Well, in case of zero frequencies, it is advised (at least I was taught so) to either merge categories, since Chi-square allows some 0 frequencies (less than 20% frequencies can be 0 when the factor has more than 2 levels), if that's reasonable, or to drop Chi-square and choose an alternative technique. As I usually drop the thing, so you may want to ask another question regarding appropriate goodness-of-fit test. $\endgroup$
    – aL3xa
    Apr 29, 2011 at 8:18
3
$\begingroup$

Would you be happier with a surgically altered goodfit object?

> 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
$\endgroup$
1
  • $\begingroup$ Cheater! /10char :) $\endgroup$ Apr 28, 2011 at 22:40
0
$\begingroup$

Try plotting it. You'll get a better idea of what's going on. As mentioned before, you're getting NaN because you're passing 0 frequencies to 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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.