Chi-squared test on small number of groups and large observations

Question

I'm comparing two data sets (3 groups) with relative large number of observations and I'm getting p-value of 0. (0 < 0.05, meaning rejecting H0, that the samples are equal, when in fact there is not a large difference)

So how I can run a sensitivity check and decide on the test sensitiveness in respect of i.e. number of observations, degrees of freedom.

I like this analysis (simulation) here provided kindly by Mr. W. Huber. But I'm not sure I can get similar sort of analysis by changing i.e. degrees of freedom or jitter the sample size etc., to obtain kind of scenario analysis.

Here is what I have tried:

Data:

 x <- data.frame(group=c(1,2,3),
            o=c(695301,154100, 224140),
            e=c(930785, 192893, 273400))


 e <- x$e  # expected frequency under Null Hypothesis
 o <- x$o # obseved frequency 
 r <- sum( (o - e) / sqrt(e))    # standardised residuals
 chq <- r^2; chq       # Chi-square statistic
 dof <- nrow(x) - 1    # degree of freedom
 p.value <- 1 - pchisq(chq, df = dof)  # p-value
 p.value  # confirm p.value with chisq.test() 

 p.val <- e/sum(e)

 chisq.test(x = o, p = p.val, rescale.p = TRUE)$p.value

Here the chisq.test has some sort of bootstrapping allowing to drop effective the degrees of freedom parameter. Here I'm getting p.value of circa 0.09.

 chisq.test(x= o, p = p.val, simulate.p.value = TRUE, B = 10)

Simple straightforward approach, yields the same result.

 chisq.test(x) # simple solution

Here I tried to perform the kind of sensitivity/or scenario analysis by simulation, without any concluding facts.

lsim <- lapply(x[ ,"o", drop=FALSE], function(x) replicate(100, jitter(x)))
o <- data.frame(lsim)
e <- x$e
apply(o, 2, function(x) { r <- sum( (x - e) / sqrt(e)); 
      chq <- sum(r^2); p.value <- 1 - pchisq(chq, df = dof); p.value } )

I would really like to see (via graphical interface) where and by which margin the two samples differ, rather than just simply concluding (without proper understanding in which group etc.) that the two sample are not equal.

Answer 1

I don't use R anywhere enough to comment on your code authoritatively, or (almost) at all, but you seem to go wrong at the outset.

Here I use Stata and trust that you can follow the results while regarding the syntax as of no interest.

Given your vectors of counts, one you are calling "observed" and the other "expected", then straight chi-square is invalid, as your vectors don't have the same total:

. chitesti  695301 154100  224140 \ 930785 192893 273400

observed frequencies from keyboard; expected frequencies from keyboard

Warning: totals of observed and expected differ
              total
observed    1.1e+06
expected    1.4e+06

The Stata command used here (for any Stata users reading this, all community-contributed commands here are from tab_chi on SSC) goes on to give output, but it should be ignored for a difference of totals that big. If what you are calling "observed" and "expected" really do have that role then you need to scale "expected" so that its sum equals the sum of "observed"; fractional expected counts are fine. Then you can repeat the chi-square test.

Alternatively, you may have a two-way table of counts, with two rows and three columns, and your misunderstanding is that one row vector should be called observed and one expected. The chi-square test is different for this set-up: rather, all counts are observed and the expected frequencies are (usually) calculated for you by your software assuming independence of rows and columns, although if there is some other hypothesis in play it determines the expected frequencies, which you need to specify.

. tabchii  695301 154100  224140 \ 930785 192893 273400

          observed frequency
          expected frequency

----------------------------------------
          |             col             
      row |        1         2         3
----------+-----------------------------
        1 |   695301    154100    224140
          | 7.07e+05  1.51e+05  2.16e+05
          | 
        2 |   930785    192893    273400
          | 9.20e+05  1.96e+05  2.81e+05
----------------------------------------

         Pearson chi2(2) = 964.1405   Pr = 0.000
likelihood-ratio chi2(2) = 962.9599   Pr = 0.000

. return list

scalars:
                  r(N) =  2470619
                  r(r) =  2
                  r(c) =  3
               r(chi2) =  964.1404788269776
                  r(p) =  4.3606891869e-210
            r(chi2_lr) =  962.959915036281
               r(p_lr) =  7.8688506561e-210

Here the chi-square test statistic is massive and the P-value minute, at most of the order of $10^{-209}$. return list here reaches into Stata's memory to pull out the maximum detail available; there is no implication from me, or Stata, that more than the first few significant figures have statistical or scientific interest.

My only other precise comment is that your R code is wrong here:

r <- sum( (o - e) / sqrt(e))    
chq <- r^2 

The problem is that chi-square is the sum of Pearson (standardized) residuals squared, not the square of the sum of Pearson residuals.

In this situation, sensitivity analysis seems moot as the result is utterly emphatic. Imagining the counts might be different makes some sense. Imagining that the number of degrees of freedom might be different makes little sense, as it is determined by the constraints of adding to a total or totals (and any other constraints). Otherwise I can't comment helpfully on the second half of the question. If the research issue is to understand why there is a big difference, looking at the residuals seems the obvious step. When P-values are minute, then any rounding by default is likely to hide important detail.

Naturally I use the only example you give, but the same principles will apply to similar examples.