2
$\begingroup$

I want to check if counts of bird territories in different cells are Poisson distributed. Due to the very small p-value produced by the below code, I have concluded they are not Poisson distributed.

Is this conclusion correct? Are the expected values and degrees of freedom in the chi-squared test correct?

(p.s. I know just by looking at a plot that they are not a Poisson distribution and I also realise the sample size might be too small to genuinely test if they are Poisson distributed. The point of this is to understand how to use the chi-squared test to check if it is Poisson distributed.)

library(dplyr);library(tidyr);library(ggplot2)

tibble(territories = c(15L, 5L, 18L, 16L, 2L, 18L, 19L, 19L, 17L, 11L, 8L, 19L, 18L, 15L, 9L, 17L, 7L, 17L, 13L, 13L, 23L, 18L, 19L, 11L, 20L, 12L, 20L, 7L, 14L, 9L, 13L, 9L, 9L, 11L, 9L, 15L, 10L)) -> df

df %>%
  dplyr::count(territories, sort = FALSE, name = "observed") ->
  df

# add zeros to data
expand.grid(territories = 2:23) %>%
  left_join(df, by = "territories") %>%
  mutate(observed = case_when(is.na(.$observed) ~ 0L, TRUE ~ .$observed)) ->
  df

df %>%
  mutate(expected = sum(df$observed) * dpois(2:23, 14)) ->
  df

stats::chisq.test(df$observed, p = df$expected, rescale.p = TRUE)

pchisq(q = 273, df= 20,lower.tail=FALSE) # manually correct degrees of freedom

df %>%
  pivot_longer(cols = observed:expected, names_to = "observed", values_to = "count") %>%
  ggplot +
  geom_col(aes(territories, count, fill = observed), position = position_dodge())

enter image description here

$\endgroup$
2
  • $\begingroup$ dplyr code is not always easy to read and the description is not very clear. What sort of counts are we looking at? What does 'territories' on the x-axis mean? $\endgroup$ Feb 6, 2023 at 11:46
  • $\begingroup$ Is territory a count variable or a label? $\endgroup$ Feb 6, 2023 at 11:48

1 Answer 1

3
$\begingroup$

The Poisson goodness-of-fit test is applied as follows. I use R as an example just as you do.

First, get the frequency table

df <- data.frame(territories = c(15L, 5L, 18L, 16L, 2L, 18L, 19L, 19L, 17L, 11L, 8L, 19L, 18L, 15L, 9L, 17L, 7L, 17L, 13L, 13L, 23L, 18L, 19L, 11L, 20L, 12L, 20L, 7L, 14L, 9L, 13L, 9L, 9L, 11L, 9L, 15L, 10L))

table(df)
territories
 2  5  7  8  9 10 11 12 13 14 15 16 17 18 19 20 23 
 1  1  2  1  5  1  3  1  3  1  3  1  3  4  4  2  1 

You see thus you have territories 2 repeated once, 5, repeated once, 7 repeated twice and so on, but there are holes in the territory values.

Fit the Poisson distribution to your data, i.e. compute the sample average:

hat_lambda <- mean(df$territories)

Now, for each territory, compute the expected counts. Each is just the probability under the estimated Poisson to observe the specified territory multiplied by the sample size. Before this, in the above table, we have to supply the other territory values that have not been observed but are probable under the estimated Poisson. There is a bit of arbitrariness in how you supply these additional values. You could bin the values in classes in order to avoid having too many classes with low frequencies.

One option could be <2, 3, 4, 6, 21, 22 and >23, which all have frequency zero and we bin the values less than 2 and those $>23$.

# sample size
N <- length(df$territories) + 7 # seven further categories
x <- sort(unique(df$territories)) # categories with positive counts
exp_counts = N*dpois(x, lambda = mean(df$territories))

Now compute the chi-squared statistic and the associated p-value

(chisq_obs <- sum((x-exp_counts)^2/exp_counts) + 7)
# p-value
1-pchisq(chisq_obs, df=length(x)+7-1)

As rightly guessed, the counts are far from being Poisson distributed, i.e. the p-value is extremely small.

Note 1. The number of degrees of freedom is equal to the number of unique observations, e.g. here the length of x, minus the number of estimated parameters, e.g. here equal to one.

Note 2. Since the extra categories provided have zero observed counts, each of them contributes to the chi-square statistic by 1. Thus a consequence of using too many categories is that you will inflate the chi-square statistic as well as the degrees of freedom of its limiting distribution.

Note 3. Since the $\chi_k^2$ is only a limiting distribution valid for large sample sizes, some writers caution the user of this approximation to be certain that the sample size is large enough s.t. each expected count is at least 5. As you can check in this example, such a rule fails. So if you want to follow this rule, you have to consider another (arbitrary) binning of the territory values in classes s.t. each class has expected value $\geq 5$.

$\endgroup$
5
  • 1
    $\begingroup$ Are zeros not needed for frequencies that are zero. For example, frequencies of 3, 4, 6 etc are missing from table(df) - are zeros not needed? $\endgroup$
    – luciano
    Feb 6, 2023 at 10:09
  • $\begingroup$ yes, you are right, we cannot have holes under the Poisson :-) I'll update my answer. $\endgroup$
    – utobi
    Feb 6, 2023 at 10:15
  • 1
    $\begingroup$ There has been much advice over several decades about minimum acceptable expected counts but experts don't agree. All expected counts at least 5 is at the hyper-cautious end of the advice. $\endgroup$
    – Nick Cox
    Feb 6, 2023 at 11:15
  • 3
    $\begingroup$ This is a careful account of the procedure, but from a scientific as well as statistical view there is much else to be said. First, Poisson may be a poor fit, but it seems unlikely that any brand-name distribution would work better for mildly irregular data such as these. Second, beyond stating that the data are not well fitted well by a Poisson the implications of that depend on what else is in mind. For example, if there are covariates, Poisson regression might still be a good idea for modelling, or at least not worse than alternatives. $\endgroup$
    – Nick Cox
    Feb 6, 2023 at 11:20
  • 1
    $\begingroup$ @NickCox, Your comments are thoughtful as usual. I'll expand my answer by adding these extra thoughts, although I'm not sure how much relevant they are to the OP's question. $\endgroup$
    – utobi
    Feb 6, 2023 at 11:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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