3
$\begingroup$

I'm confused about the degrees-of-freedom discrepancy between R's chisq.test function, and doing the calculation "by hand".

Using the code below, the test statistics in the two approaches agree exactly but chisq.test reports 11 degrees-of-freedom, whereas doing it "by hand" I would have used k-p-1=10 where k=number of categories (12), p = number of parameters estimated from the data (1 on account of estimating the mean). I don't see how the chisq.test can know the number of parameters estimated from the data and how the number of degrees-of-freedom can agree?

k <- 0:11
obs <- c(57,203,383,525,532,408,273,139,45,27,10,6)

#Get mean...
x <- rep(k,times=obs)
lambda <- mean(x)

#probabilites under H0 (Poisson model) and expected values ek...
p <- dpois(0:10,lambda)
# p(k >= 11)
prest <- 1-sum(p)
pall <- c(p,prest)
ek <- pall*sum(obs)

#R test...
xsq <- chisq.test(obs,p=pall)

#Manual...
xsqts <- sum((obs-ek)**2/ek)
xsqdf <- length(k)-1-1
pval <- 1-pchisq(xsqts,xsqdf)
```
$\endgroup$
7
  • 2
    $\begingroup$ (1) The first line of your code is a bomb waiting to explode in the face of the unwary. Kindly remove it. (2) There are multiple errors in both calculations, beginning with setting the expected number of values exceeding $10$ to zero (it's almost $6$). $\endgroup$
    – whuber
    Aug 10 at 21:01
  • 2
    $\begingroup$ I second the removal of rm(list=ls()). $\endgroup$
    – Dave
    Aug 10 at 21:05
  • 1
    $\begingroup$ rm(list=ls()) removed. Why is the expected number of values exceeding 10 zero? Here are the values "ek". Doesn't this show it's 5.77? [1] 54.418655 210.580164 407.433861 525.539688 508.411286 393.472906 [7] 253.765885 140.282938 67.855416 29.175054 11.289672 5.774475 $\endgroup$
    – PM.
    Aug 10 at 21:24
  • 1
    $\begingroup$ The point is there's no point at which you tell it you did that (and no argument to chisq.test where you could have) - it obviously isn't written to cover that situation. You could use it to calculate the statistic and then look up pchisq yourself, just as you did, but if it's something you need often, you could write a function that does it. The big problem is that if you bin together values at all (as you must at least do at the upper end with the Poisson), your parameter estimate has to be based on the binned (censored) values, not the original ones, or the statistic isn't chisquared. $\endgroup$
    – Glen_b
    Aug 11 at 11:01
  • 1
    $\begingroup$ ... even asymptotically. Then there's the second problem that the choice of bins if you do bin shouldn't be based on the observed counts, which is difficult to arrange. $\endgroup$
    – Glen_b
    Aug 11 at 11:02
2
$\begingroup$

Besides the points highlighted in the comments, I guess the issue is with the way degrees of freedom is calculated in chisq.test when x (here, obs) is a vector. There is an explanation in the documentation:

If x is a matrix with one row or column, or if x is a vector and y is not given, then a goodness-of-fit test is performed (x is treated as a one-dimensional contingency table). The entries of x must be non-negative integers. In this case, the hypothesis tested is whether the population probabilities equal those in p, or are all equal if p is not given.

This wikipage explains how to calculate the degrees of freedom for a goodness-of-fit test:

For a test of goodness-of-fit, df = Cats − Parms, where Cats is the number of observation categories recognized by the model, and Parms is the number of parameters in the model adjusted to make the model best fit the observations: The number of categories reduced by the number of fitted parameters in the distribution.

So, in this case, it should be length(k)-1, 11. And this is the number reported by chisq.test.

Edit:

Snedecor and Cochran (1989, p.77) explains the calculation:

The theoretical frequency distribution of $s^2$ under normality depends on the population variance $\sigma^2$ and on the degrees of freedom $\nu$ in $s^2$. Since both of these were known when the theoretical distribution in table 5.9.1 [a table with frequencies and information about $\sigma^2$] was fitted, we say that the theoretical distribution is completely specified. In this case, the number of degrees of freedom in $\chi^2$ is $(k — 1)$, where $k$ is the number of classes whose contributions are summed in finding $\chi^2$- Note that the degrees of freedom in $\chi^2$ do not involve the sample size $n$ [emphasis mine].

In the next paragraph, they continue:

If we want to test whether this theoretical distribution fits but do not know the value of $\sigma^2$, we can use the average of the 511 sample variances $s^2$ [in the above mentioned table] as an estimate of $\sigma^2$ when fitting. The number of degrees of freedom in $\chi^2$ is then $(k —2)$…. Similarly, in fitting a normal distribution, the degrees of freedom in $\chi^2$ are $(k — 1)$ if $\mu$ and $\sigma$ are both known for the populations when fitting,…,but the degrees of freedom are $(k — 3)$ if $\bar X$ and $s$ from the sample are used as estimates of $\mu$ and $\sigma$ when fitting. The rule is

$$\text{df = (k - 1 - number of fitted parameters)}$$

So, answer to your last question (as @Glen_b pointed out): it can’t. For example, if you use pearson.test from nortest package you can see that it calculates ($k-3$) degrees of freedom (if you specify n.classes and adjust=TRUE) because sample is used to estimate parameters (but there is a warning, see the note). This is not the case for chisq.test. [I hope these long quotes don't make the answer difficult to read.]

$\endgroup$
2
  • $\begingroup$ I'm confused as that doesn't seem to match this which says dof = k-c where c = number of params estimated + 1. There's a +1 there? itl.nist.gov/div898/handbook/eda/section3/eda35f.htm $\endgroup$
    – PM.
    Aug 11 at 8:13
  • 1
    $\begingroup$ I checked the Snedecor and Cochran (1989) reference (in the link) and I think I understand the source of confusion. It is about having theoretical distribution completely specified (or not). I will extend my answer when I have time. $\endgroup$
    – T.E.G.
    Aug 11 at 9:04

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.