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:


 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.

  • $\begingroup$ "in fact there is not a large difference": why do you say this? Absolutely, the counts differ a lot. Relatively, which is what the chi-square test assesses on one scale, the rows or columns differ a lot. What scale or standard are you using there? $\endgroup$ – Nick Cox Oct 4 '18 at 18:16
  • $\begingroup$ I guess I have not picked the right example, this group is not that representative. I have cases where the sample differs relatively slightly and I get p.values close to zero. $\endgroup$ – Maximilian Oct 4 '18 at 19:15
  • $\begingroup$ Can please someone at least confirm that the approach/solution above is correct? Are there any errors? Many thanks $\endgroup$ – Maximilian Oct 5 '18 at 6:45
  • $\begingroup$ There is a meta-issue of why your question got so little attention, as chi-square tests for categorical counts are not esoteric and people fluent in R the biggest active group of those who answer questions. I have two guesses that may help with your future questions. Sometimes bundling several related questions together makes a question too complicated to appeal or to be easy to answer, especially if there appears to be some confusion about principles. Also, any flavour of "please try to follow my code and comment on it" often deters, even if what you are doing statistically is the issue. $\endgroup$ – Nick Cox Oct 5 '18 at 8:06
  • $\begingroup$ Thanks so much Nick, I really appreciate it, also your answer (+1). Yes, there is indeed bug in my code, with the summation. If you look into the simulation function, I use correctly summation sum(r^2) but actually use the summation twice which is also an error. I will read your answer correctly and I'm hopeful to guide me further. Many thanks. $\endgroup$ – Maximilian Oct 5 '18 at 8:08

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
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

                  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.


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.