7
$\begingroup$

We have a big urn that contain $N_{Tot}$ balls. Balls are of $r$ different colours. The number of balls of the $i^{th}$ colour (before sampling) is $N_i$. John sampled $x$ balls in total (without replacement) from this urn. The number of balls of each colour are $x_1, x_2, ..., x_{r-1}, x_r$.

The question I am trying to answer: Did John sample balls independently of the number of balls existing in it's colour?


My thoughts

I first thought a simple regression with $r$ data points, where the x variable is the number of balls in this color and the y coordinate is the number of balls that were sampled in this color. But that wouldn't work. One obvious reason is that each point is not an independent observation as sampling the first ball from one colour means not sampling the first ball in another colour.

$x_i$ should follow a hypergeometric distribution with parameters $N_{tot}-N_i$ (failure), $N_i$ (success) and $x$ (total number of balls drawn), where $N_{tot} = \sum_{i=1}^r N_i$ is the total number of balls. I can compare $x_i$ to the median of the hypergeometric distribution to tell for each colour whether the urn was rather overrepresented in the sample or underrepresented. Let's call this binary variable $B$

Then, I can perform a Spearman correlation between $B$ and the number of balls in each colour $N_1, N_2, ..., N_{r-1}, N_r$. If the p.value of the Spearman correlation is low enough then, I can conclude that the balls were not drawn at random independently of the colour of the balls.

Does it seem to be a good methodology?

Would it be a test with low power or would it be completely wrong?

$\endgroup$
7
  • 2
    $\begingroup$ +1 This is an example of the multivariate hypergeometric distribution. $\endgroup$
    – whuber
    Jun 15, 2015 at 22:21
  • $\begingroup$ Thnks that seems interesting. I don't quite get how I can ask whether the number of balls in the colour changes the probability of sampling one given ball of that color (Ha) or not (Ho). $\endgroup$
    – Sulawesi
    Jun 15, 2015 at 23:12
  • $\begingroup$ I would love to give a bounty to attract more attention but I am afraid I don't have enough reputation for that! $\endgroup$
    – Sulawesi
    Jun 19, 2015 at 17:47
  • 1
    $\begingroup$ Is there some reason why you can't use a Chi-square test for this? It seems that you can calculate the expected number of each color of ball collected, based on the multivariate hypergeometric distribution. Then a standard Chi-square test based on observed numbers of each color would indicate whether the sample was significantly far from what would have been expected. $\endgroup$
    – EdM
    Jun 19, 2015 at 19:01
  • 1
    $\begingroup$ Ok. I think you can make an answer out of your comments. I am still a bit confused about the non-independence of the data points. If there are two colors only and John draws 10 balls in total, including 8 yellow balls, then he cannot draw more than 2 blue balls. in other words, the $\chi^2 = \frac{(8 - \bar x_1)^2}{\bar x_1} + \frac{(3 - \bar x_2)^2}{\bar x_2}$ is not a possible outcome. I feel like this inter-dependence would not allow one to perform a Chi-square test. $\endgroup$
    – Sulawesi
    Jun 25, 2015 at 0:04

1 Answer 1

4
$\begingroup$

Because your proposed test has no theoretical foundation and does not account for the correlations among the counts, it would not be a good use of our time to evaluate it. Instead, let's develop some tests that will work.

The Likelihood Ratio test is guaranteed to work well when all counts are relatively large. So will the chi-squared test provided only a small proportion of all balls in the urn are taken in the sample. (The usual rule of thumb is to be cautious when the sample exceeds ten percent of the total.) In other cases, simulating from the null distribution is effective.

The likelihood of observing counts $x_1,\ldots, x_r$ in a sample from an urn with $N_1,\ldots, N_r$ balls of each color is

$$\mathcal{L}(x;N) = \binom{N_1}{x_1}\binom{N_2}{x_2}\cdots \binom{N_r}{x_r}\ /\ \binom{N_1+\cdots+N_r}{x_1+\cdots+N_r}.$$

This is best expressed in terms of its deviance $D=-2\log\mathcal L$ because asymptotically, the deviance has a $\chi^2$ distribution with $r-1$ degrees of freedom

By sampling a few thousand times (with the computer) you can estimate the distribution of $D.$ The p-value of the observed values is the area of the right tail determined by the deviance of the observations.

As an example, consider an urn with $r=7$ colors in the quantities $(N_i) = (34,45,41,35,49,47,51,42).$ I computed the null deviance in ten thousand samples of size $167$ (half the balls are in each sample), shown in the figure at the left. I also computed the deviance in another ten thousand samples where the colors were selected with different probabilities ranging from 20% (for color $1$) down to just 6% (for color $7$). The latter deviances tend to be large. Indeed, approximately 96% of them exceed the 95th percentile of the null deviances. In this sense, the power of this test, when conducted at the $alpha=100-95\% = 5\%$ level, is 96%.

Figure 1

We can play the same game with the usual $\chi^2$ test. If this test were directly applicable, the null distribution of its p-values would be uniform. The null distribution is skewed, though, with most of the p-values very large: see the left graphic in the next plot. But, upon examining this simulation, we could declare a chi-squared test result to be significant at the $\alpha$ level whenever its p-value is less than the $\alpha$ quantile of this null distribution. That quantile is near $0.4$ (rather than the expected $0.05$), shown with the shaded bars.

The right hand graphic similarly displays the distribution of chi-squared test p-values under the alternative hypothesis. Its "nominal power" is the rate at which the reported (nominal) p-value is less than $\alpha,$ shown with the area of the red bar. That is only 48.2%, far less than the 96% achieved by the Likelihood Ratio (LR) test. However, when we use the $\alpha$ quantile of the null distribution as our threshold, now the null hypothesis is rejected in 95% of the simulated data. The power, 95%, is essentially the same as the LR test.

Figure 2

What these simulations demonstrate, then, is

  1. It is invalid to apply the standard chi-squared test or LR test unless the sample is a small portion of the urn and the sample counts are fairly large. The skewed null distribution in the second figure shows what goes wrong.

  2. Nevertheless, both of these tests can be used provided the p-value is computed using the actual null distribution (as estimated through simulation) rather than using the standard formulas (which rely on asymptotically large samples of urns with huge populations).

The R code to perform these computations and display these figures follows.

#
# Multivariate hypergeometric distribution.
# Draw `n` balls without replacement from an urn with length(N) colors, each
# appearing N[i] times.
#
# Returns a vector of counts of the colors.
#
rmhyper <- function(n, N, p) {
  if (missing(p)) p <- rep(1/sum(N), length(N))
  prob <- rep(p, N)
  prob <- prob / sum(prob)
  x <- sample(rep(seq_along(N), N), n, prob=prob)
  tabulate(x, length(N))
}
#
# Returns the likelihood of observing `k` in a hypergeometric draw specified
# by `N`.
#
dmhyper <- function(k, N, log.p=TRUE) {
  q <- sum(lchoose(N, k)) - lchoose(sum(N), sum(k))
  if (!isTRUE(log.p)) q <- exp(q)
  return(q)
}
#
# Perform a chi-squared test.
#
mhyper.test <- function(k, N, ...) {
  chisq.test(k, p=N / sum(N), ...)
}
#
# Simulations.
#
alpha <- 0.05
n.sim <- 1e4
set.seed(17)
N <- rpois(8, 40) + 1          # Determine the population randomly
p <- rgamma(length(N), 10/4)   # Determine the alternative hypothesis randomly
p <- rev(sort(p / sum(p)))
n <- ceiling(sum(N)*0.5)
#
# Display the sampling probabilities.
#
plot(p, type="h", col=rainbow(length(N), .9, .8),
     lwd=3, ylim=c(0, max(p)), xlab="Color", main="Probabilities")
#
# The chi-squared test p-values are ok when the sample is a small fraction
# of the population.
#
p.values.null <- replicate(n.sim, mhyper.test(rmhyper(n,N), N, simulate.p.value=FALSE)$p.value)
p.values.alt <- replicate(n.sim, mhyper.test(rmhyper(n,N, p), N, simulate.p.value=FALSE)$p.value)
power <- round(100*mean(p.values.alt <= alpha), 1)
power.alt <- round(100*mean(p.values.alt <= quantile(p.values.null, alpha)))
k <- round(20 * quantile(p.values.null, alpha)) - 1

b <- seq(0, 1, by=0.05)
par(mfrow=c(1,2))
h <- hist(p.values.null, freq=FALSE, breaks=b, xlab="p",
     col=c("#d08080", rep("#e0e0e0", k), rep("White", 20-k-1)),
     main="Histogram of Chi-square P-values\nNull Distribution")
hist(p.values.alt, freq=FALSE, breaks=b, xlab="p", ylim=c(0, max(h$density)),
     col=c("#d08080", rep("#e0e0e0", k), rep("White", 20-k-1)),
     sub=bquote(paste("Nominal power = ", .(power), "%; Simulation power = ", .(power.alt), "%")),
     main="Histogram of Chi-square P-values\nAlternative Distribution")

q.null <- -2 * replicate(n.sim, dmhyper(rmhyper(n, N), N))
q <- -2 * replicate(n.sim, dmhyper(rmhyper(n, N, p), N))


xlim <- range(c(q, q.null))

h <- hist(q.null, xlim=xlim, freq=FALSE, breaks=50, xlab=expression(-2~~log(p)),
     col="#f0f0f0",
     main="Null distribution of deviance")
abline(v = quantile(q.null, 1-alpha), col="Red", lwd=2)

power <- signif(mean(q >= quantile(q.null, 1-0.05)), 2)
hist(q, xlim=xlim, freq=FALSE, breaks=50, xlab=expression(-2~~log(p)),
     col="#f0f0f0", ylim=c(0, max(h$density)),
     main="Alternative distribution of deviance",
     sub=bquote(paste("Power = ", .(100*power), "%")))
abline(v = quantile(q.null, 1-alpha), col="Red", lwd=2)
par(mfrow=c(1,1))
$\endgroup$

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.