5
$\begingroup$

Let's say I have 7 urns filled with random numbers of colourful marbles. An example dataset is as below:

data <- matrix (c(5,3,4,4,4,2,1,1,1,1,1,2,2,1,2,2,1,1,2,1,4,1,1,2,4,1,3,1,7,1,1,2,1,3,3), ncol = 5);
rownames(data) <- as.character(seq(1,7));
colnames(data) <- c("red", "blue", "yellow", "green", "pink");

I assume the count-based colour distribution to be multinomial. I want to test whether the contents of Urn 1 follow the colour distribution of Urns 2-7.

I get the MLEs for Urns 2-7 as:

p_sample <- colSums(data[2:7,])/ sum(colSums(data[2:7,]))

1) When I run a $\chi^2$ test as below, R displays a warning message since the expected counts (EC)<5.

obs <- data["1",]
chisq.test(obs,p=p_sample)
# Chi-squared test for given probabilities

# data:  obs
# X-squared = 8.0578, df = 4, p-value = 0.08948
#
# Warning message:
# In chisq.test(obs, p = p_sample) :
#  Chi-squared approximation may be incorrect

One of the answers to this question states that $\chi^2$ test would nevertheless return accurate results as long as ECs exceed 1.0 if a very simple $\frac{N-1}{N}$ correction is applied to the test statistic. Is the correction implemented simply like this: $\chi^2 = (\sum_i \frac{(O_i - E_i)^2}{E_i} )* \frac{N-1}{N} $ ?

This link suggests so but I am not sure about its reliability:

If one has the regular Pearson chi-square (e.g., in the output from statistical software), it can be converted to the 'N - 1' chi-square as follows:

           'N -1' chi-square = Pearson chi-square x (N -1) / N

2) Alternative to $\chi^2$, since the number of counts is low, I ran a Fisher's exact test and got the results below. Is how I call the fisher.test below correct? The result I get makes me think: "No". I am confused since the R help document only refers to use cases for contingency tables.

fisher.test(obs, sum(obs)*p_sample)

# Fisher's Exact Test for Count Data
#
# data:  obs and sum(obs) * p_sample
# p-value = 1
# alternative hypothesis: two.sided
$\endgroup$
  • $\begingroup$ You should't use the same observation twice. You probably want to compare the color distribution in urn 1 with the color distribution of the other urns combined. The null hypothesis of the corresponding (generalized) Fisher's exact test is 'the two distributions are equal' which is similar to 'the distribution in urn 1 is equal to the one of all urns together'. $\endgroup$ – Michael M Nov 15 '13 at 17:01
  • $\begingroup$ @MichaelMayer, thanks I will update the code and question accordingly to get answers to 2) and 3) $\endgroup$ – Zhubarb Nov 15 '13 at 17:04
6
$\begingroup$

As already pointed out in my comment referring to the original question, your preferred null hypothesis "color distribution in Urn 1 is equal to color distribution in all Urns combined" is equivalent to the null hypothesis "color distribution in Urn 1 is equal to color distribution in Urn 2-7". The former recycles observations in Urn 1, destroying independence. Thus, for technical reasons, the latter formulation is better suited.

Null hypotheses like this are usually tested by the approximate $\chi^2$-test of independence and, for not too large samples, by Fisher's exact test for general two-dimensional tables. Your sample is not large, so in principle we can go with Fisher's version and do not need to care about improving approximation by $\chi^2$. Nevertheless, we can compare the different ways in your example and see what comes up:

Fisher's exact test for general tables:

data2 <- cbind(data[1,], colSums(data[2:7,]))
data2

# Data looks like this
       [,1] [,2]
red       5   18
blue      1    8
yellow    2   11
green     1   12
pink      7   11

fisher.test(data2)  # Yields: p value 0.2977

Quite the same with the classic $\chi^2$-test:

chisq.test(data2)

# Output
X-squared = 5.6345, df = 4, p-value = 0.2282

Warning message:
In chisq.test(data2) : Chi-squared approximation may be incorrect

Now, the slightly less optimistic version with $(N-1)/N$ correction:

N <- sum(data2)
1 - pchisq(5.6345*(N-1)/N, 4) # Yields a p value of 0.234

Or, calculating degrees of freedom and Pearson's statistic from scratch

N <- sum(data2)                     # Number of obsevations
E <- outer(rowSums(data2), colSums(data2), "*")/N   # Expected cell counts
T <- sum((data2-E)^2/E)                 # Pearson's statistic
df <- prod(dim(data2)-1)                # degrees of freedom
1 - pchisq(T*(N-1)/N, df)                           # Corrected p value 

Thus, in your example (no matter what procedure we use) and at the 5% level, we cannot claim that color distribution in Urn 1 is different to the rest. We also cannot claim that the color distribution in Urn 1 is equivalent to the rest.

The discussion and references in the link in your post show advantages and disadvantages of the three procedures (among others). Given that the procedure needs to be simple, it seems to be reasonable to use the $N-1$-corrected version of the original $\chi^2$-test of independence. Please make sure to use the formula for two-dimensional tables, i.e. the test statistic of the $\chi^2$-test of independence (http://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test).

What went wrong in your two attempts?

$\chi^2$: You performed a $\chi^2$ goodness of fit test. This is wrong as you do not know the true color distribution of Urns 2-7. You only know the empiric distribution.

Fisher's exact test: You multiply the relative counts in Urns 2-7 with the sample size of Urn 1, which doesn't make much sense.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you Michael, this is very helpful. So in this case can we not apply a $\frac{N-1}{N}$ correction to ensure the $\chi^2$ test produces accurate results? I am investigating this as part of a model that is supposed to run on a web script - which may not be able to access the FortRan and C code that runs behind the fisher.test function in R. If it were an option the $\chi^2$ test would be more convenient implementation-wise. $\endgroup$ – Zhubarb Nov 16 '13 at 11:42
  • 1
    $\begingroup$ Modified the answer regarding this point. I wasn't aware that computational simplicity is required. $\endgroup$ – Michael M Nov 16 '13 at 12:46

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.