Compare chi-squared goodness of fit results for Poisson and Negative binomial

Question

Steps followed for testing :

1. Derive parameter estimates using fitdistr() function in MASS package for the dataset, dt consisting of discrete values varying from 0 to 50.

df=dt[[1]]
nbfit <- fitdistr(df,'negative binomial')
pfit <- fitdist(df, "pois")

2. Use parameter estimates for creating reference distribution

    fitnb <- dnbinom(0:50, size=0.4, mu=3.9)
    fitp <- dpois(0:50, lambda=3.9)

I found that lambda value is same as mu.

3. Get frequencies

    t <- table(df)   
    D <- as.data.frame(t)
    observed_freq <- D$Freq

4. perform the chi-squared test

chisq.test(observed_freq, fitnb, simulate.p.value = TRUE)
chisq.test(observed_freq, fitp , simulate.p.value = TRUE)

5. Results

for NB

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  observed_freq and fitnb
X-squared = 476, df = NA, p-value = 0.989

for Poisson

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  observed_freq and fitp
X-squared = 476, df = NA, p-value = 0.9875

Question:

1. What can we comment upon which distribution fits better? Clearly, for both X-squared is the same and we fail to reject the null hypothesis.

2. Can I perform some other test?

Edit

I also used lrtest() function and got the following result:

fit_poi <- fitdistr(df,"poisson")
fit_nbin <- fitdistr(df,"negative binomial")
lrtest(fit_poi,fit_nbin)

Result

Model 1: fit_poi
Model 2: fit_nbin
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   1 -2537.7                         
2   2 -1159.3  1 2756.6  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

With this, can we say negative binomial performs better?

Histogram of actual data used in the test (Note: Codes differ a little based on this dataset)

Actual data

value   count
0   205
1   65
2   40
3   40
4   32
5   18
6   15
7   11
8   8
9   9
10  7
11  6
12  8
13  3
14  1
15  8
16  2
18  2
19  2
20  5
21  1
23  2
24  1
25  1
29  1
31  1
32  2
35  1
36  1
42  1
43  1
45  1
50  1
53  1

Answer 1

About 1. Your chi-squared test does not reject $H_0$, thus the two fits are statistically equivalent; By Occam's Razor principle the Poisson model, i.e. the most parsimonious of the two, is to be preferred. However, as noted in the comments by whuber, these tests seem to give strange results; see the Addendum below.

About 2. The test you are looking for is sometimes called a test for overdispersion since the Poisson distribution imposes the "mean equals variance" assumption. Since the Poisson distribution is a particular case of the Negative Binomial distribution when the latter is expressed in a particular parametrization, then a likelihood ratio test can be applied. However, this test doesn't have the usual $\chi^2$ limiting distribution due to the fact that the parameter being tested is assumed on the border of the parameter space.

The easiest solution is then to apply a score test (Dean and Lawless, 1989 "Tests for detecting overdispersion in Poisson regression models", Journal of the American Statistical Association 84: 467–472.). This amounts to applying a $t$-test to

$$ Z_i = \frac{(Y_i - \mu_i)^2 - \mu_i}{\mu_i\sqrt{2}}. $$

The test is post-hoc, in the sense that it is performed subsequent to modelling the data.

Here is an R example applied to a simulated variable from the Poisson distribution.

# generate some data
set.seed(12)
x <- data.frame(y = rpois(30,3.5))

fit_poi <- glm(y~1, family = "poisson", data = x)
mu <- predict(fit_poi, type="response")
z <- ((x$y - mu)**2 - x$y)/ (mu * sqrt(2))
zscore <- lm(z ~ 1)
summary(zscore)

Call:
lm(formula = z ~ 1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.67076 -0.55149 -0.27887  0.07895  2.87330 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.08207    0.16413    -0.5    0.621

Residual standard error: 0.899 on 29 degrees of freedom

The p-value is $0.621$ thus we cannot reject the null hypothesis and conclude that the data are not overdispersed.

Addendum

With your data the score test gives:

summary(zscore)
Call:
lm(formula = z ~ 1)

Residuals:
   Min     1Q Median     3Q    Max 
 -9.19  -8.19  -5.72  -5.72 422.01 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    8.463      1.660   5.099 4.85e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 37.22 on 502 degrees of freedom

The conclusion is that the data are definitely overdispersed, so non-Poisson. You should then opt for the Negative Binomial and as suspected above (and by whuber in the comments to your post) there must be something wrong with your first chi-squared test.