Determine if three is statistically different than ten for a very large number of observations (1,000,000)

Question

For 1,000,000 observations, I observed a discrete event, X, 3 times for the control group and 10 times for the test group. How do I determine for a large number of observations (1,000,000), if three is statistically different than ten?

Answer 1

I think a simple chi-squared test will do the trick. Do you have 1,000,000 observations for both control and test? If so, your table of observations will be (in R code)

Edit: Woops! Left off a zero!

m <- rbind(c(3, 1000000-3), c(10, 1000000-10))
#      [,1]   [,2] 
# [1,]    3 999997
# [2,]   10 999990

And chi-squared test will be

chisq.test(m)

Which returns chi-squared = 2.7692, df = 1, p-value = 0.0961, which is not statistically significant at the p < 0.05 level. I'd be surprised if these could be clinically significant anyway.

Answer 2

The huge denominators throw off one's intuition. Since the sample sizes are identical, and the proportions low, the problem can be recast: 13 events occurred, and were expected (by null hypothesis) to occur equally in both groups. In fact the split was 3 in one group and 10 in the other. How rare is that? The binomial test answers.

Enter this line into R: binom.test(3,13,0.5,alternative="two.sided")

The two-tail P value is 0.09229, identical to four digits to the results of Fisher's test.

Looked at that way, the results are not surprising. The problem is equivalent to this one: If you flipped a coin 13 times, how surprising would it be to see three or fewer, or ten or more, heads. One of those outcomes would occur 9.23% of the time.

Answer 3

A (two-sided) Fisher's Exact test gives p-value = 0.092284.

function p = fexact(k, x, m, n)
%FEXACT Fisher's Exact test.
%   Y = FEXACT(K, X, M, N) calculates the P-value for Fisher's
%   Exact Test.
%   K, X, M and N must be nonnegative integer vectors of the same
%   length.  The following must also hold:
%   X <= N <= M, X <= K <= M and K + N - M <= X.  Here:
%   K is the number of items in the group,
%   X is the number of items in the group with the feature,
%   M is the total number of items,
%   N is the total number of items with the feature,

if nargin < 4
   help(mfilename);
   return;
end
nr = length(k);
if nr ~= length(x) | nr ~= length(m) | nr ~= length(n)
   help(mfilename);
   return;
end

na = nan;
v = na(ones(nr, 1));
mi = max(0, k + n - m);
ma = min(k, n);

d = hygepdf(x, m, k, n) * (1 + 5.8e-11);
for i = 1:nr
  y = hygepdf(mi(i):ma(i), m(i), k(i), n(i));
  v(i) = sum(y(y <= d(i)));
end
p = max(min(v, 1), 0);
p(isnan(v)) = nan;

For your example, try fexact(1e6, 3, 2e6, 13).

Answer 4

In this case Poisson is good approximation for distribution for number of cases. There is simple formula to approximate variance of log RR (delta method) .

log RR = 10/3 = 1.2, se log RR = sqrt(1/3+1/10) = 0.66, so 95%CI = (-0.09; 2.5)

It is not significant difference at 0.05 level using two-sided test. LR based Chi-square test for Poisson model gives p=0.046 and Wald test p=0.067. This results are similar to Pearson Chi-square test without continuity correction (Chi2 with correction p=0.096). Another possibility is chisq.test with option simulate.p.value=T, in this case p=0.092 (for 100 000 simulations).

In this case test statistics is rather discrete, so Fisher test can be conservative. There is some evidence that difference can be significant. Before final conclusion data collecting process should be taken into account.

Answer 5

If you wanted to check non-parametrically for significance, you could bootstrap the confidence intervals on the ratio, or you could do a permutation test on the two classes. For example, to do the bootstrap, create two arrays: one with 3 ones and 999,997 zeros, and one with 10 ones and 999,990 zeros. Then draw with replacement a sample of 1m items from the first population and a sample of 1m items from the second population. The ratio we're interested in is the ratio of "hits" in the first group to the ratio of "hits" in the second group, or: (proportion of ones in the first sample) / (proportion of ones in the second sample). We do this 1000 times. I don't have matlab handy but here's the R code to do it:

# generate the test data to sample from
v1 <- c(rep(1,3),rep(0,999997))
v2 <- c(rep(1,10),rep(0,999990))
# set up the vectors that will hold our proportions
t1 <- vector()
t2 <- vector()
# loop 1000 times each time sample with replacement from the test data and
# record the proportion of 1's from each sample
# note: this step takes a few minutes. There are ways to write it such that
# it will go faster in R (applies), but it's more obvious what's going on this way:
for(i in 1:1000) {
   t1[i] <- length(which(sample(v1,1000000,replace=TRUE)==1)) / 1000000
   t2[i] <- length(which(sample(v2,1000000,replace=TRUE)==1)) / 1000000
}
# what was the ratio of the proportion of 1's between each group for each random draw?
ratios <- t1 / t2
# grab the 95% confidence interval over the bootstrapped samples
quantile(ratios,c(.05,.95))
# and the 99% confidence interval
quantile(ratios,c(.01,.99))

The output is: 5% 95% 0.0000000 0.8333333 and: 1% 99% 0.00 1.25 Since the 95% confidence interval doesn't overlap the null hypothesis (1), but the 99% confidence interval does, I believe that it would be correct to say that this is significant at an alpha of .05 but not at .01.

Another way to look at it is with a permutation test to estimate the distribution of ratios given the null hypothesis. In this case you'd mix the two samples together and randomly divide it into two 1,000,000 item groups. Then you'd see what the distribution of ratios under the null hypothesis looks like, and your empirical p-value is how extreme the true ratio is given this distribution of null ratios. Again, the R code:

# generate the test data to sample from
v1 <- c(rep(1,3),rep(0,999997))
v2 <- c(rep(1,10),rep(0,999990))
v3 <- c(v1,v2)

# vectors to hold the null hypothesis ratios
t1 <- vector()
t2 <- vector()

# loop 1000 times; each time randomly divide the samples
# into 2 groups and see what those two random groups' proportions are
for(i in 1:1000) {
  idxs <- sample(1:2000000,1000000,replace=FALSE)
  s1 <- v3[idxs]
  s2 <- v3[-idxs]
  t1[i] <- length(which(s1==1)) / 1000000
  t2[i] <- length(which(s2==1)) / 1000000
}

# vector of the ratios
ratios <- t1 / t2

# take a look at the distribution
plot(density(ratios))

# calculate the sampled ratio of proportions
sample.ratio <- ((3/1000000)/(10/1000000))
# where does this fall on the distribution of null proportions?
plot(abline(v=sample.ratio))

# this ratio (r+1)/(n+1) gives the p-value of the true sample
(length(which(ratios <= sample.ratio)) + 1) / (1001)

The output is ~ .0412 (of course this will vary run to run since it's based on random draws). So again, you could potentially call this significant at the .05 value.

I should issue the caveats: it depends too on how your data was collected and the type of study, and I'm just a grad student so don't take my word as gold. If anyone has any criticism of my methods I'd love to hear them since I'm doing this stuff for my own work as well and I'd love to find out the methods are flawed here rather than in peer review. For more stuff like this check out Efron & Tibshirani 1993, or chapter 14 of Introduction to the Practice of Statistics by David Moore (a good general textbook for practitioners).

Answer 6

I would be really surprised if you find the difference statistically significant. Having said that you may want to use a test for a difference of proportions (3 out of 1M vs 10 out of 1M).

Answer 7

In addition to the other answers: If you have 1,000,000 observations and when your event comes up only a few times, you are likely to want to look at a lot of different events. If you look at 100 different events you will run into problems if you work with p<0.05 as criteria for significance.