1
$\begingroup$

What is the p-value for each individual test in my prop.test function? (See code below).

When doing multiple testing (k = 10 000 tests in this case), I want to find the alpha for each individual test in prop.test, since it obviously is not 0.05. The function prop.test must be adjusting the p-value somehow.

Upon reading on the internet and in R (typing prop.test in the console), I did not find an answer to what the p-value is and how to obtain it. I realize that typing prop.test in the console can help but I tried it and could not understand the prop.test-code well enough to understand what the alpha (for individual test) is. Hence, I would appreciate if someone could explain how I can obtain this number for alpha (for individual tests).

Thanks a lot in advance, /Pedram

(CODE:)

k <- 10000

pH0 <- 1:k
pH1 <- 1:k
nA <- 4000
nB <- 4000
p0 <- 0.01
p1 <- 2*p0

a = 0.05
pD = a

for(i in 1:k)
{
      x1 <- rbinom(1,nA,p0)
      x2 <- rbinom(1,nB,p0)
      y <- rbinom(1,nA,p1)

      pH0[i] <- prop.test(c(x1,x2),c(nA,nB))$p.value < pD
      pH1[i] <-prop.test(c(x1,y),c(nA,nA))$p.value < pD 
}`
$\endgroup$
  • $\begingroup$ Can you please elaborate what the problem is? You can easily set the level for the confidence interval by adding conf.level = 0.95. The p-value, however, is independent from the choice of alpha, because (as always) it is simply the probability of your data under the null hypothesis. In your code, prop.test has no chance to correct for multiple testing, because you are simply repeating the test 10000 times without telling R that you want to do multiple testing. $\endgroup$ – LuckyPal Apr 16 '19 at 8:03
  • $\begingroup$ Thank you for your comment. I will try to elaborate. I am simulating 10 000 tests, I want to simulate 10 000 "rows" of data for 4000 healthy individuals (these individuals are my first binomial sample) and for 4000 unhealthy individuals (my second binomial sample). I want to test this and to find try and find a cut-off value (for example, if p-value < 0.0005) then we find significance to make my FDR/FWER = alpha. Is it clearer now? $\endgroup$ – Pedram Rayat Apr 16 '19 at 8:08
  • $\begingroup$ Thank you, now it's more clear to me what you are doing! $\endgroup$ – LuckyPal Apr 16 '19 at 10:52
0
$\begingroup$

If I understood correctly what you are doing, the problem is basically that sum(pH0)/k is not 0.05, is that right?

This problem can easily be solved: just add the option , correct = F to prop.test, i.e. in your for-loop you should write: pH0[i] <- prop.test(c(x1,x2),c(nA,nB), correct = F)$p.value < pD

By default, prop.test uses Yate's correction for continuity, which might be appropriate in situations where sample sizes are small (or to be more precise, when the expected frequencies are small). However, this correction tends to overcorrect and leads to type-I-error rate smaller than desired, and is not appropriate for samples as large as yours.

Edit based on OP's comment

Maybe now I got what you are trying to achieve. The issue of multiple testing is not simply performing multiple tests. Actually, the whole idea of p-values in the frequentist framework relies on multiple tests, i.e. the idea that you can potentially repeat the experiment infinitely many times.

Instead, multiple testing becomes a problem if one single decision is based on multiple hypothesis tests. The following code is quite inelegant but I hope it is at the same time very explicit to you in what it is doing.

At the bottom you see that each single sub-hypothesis test fulfills the 5% requirement. However, if we decide to reject the null hypothesis whenever one of the sub-hypotheses was rejected, re have increased type I error rate. To be precise: with three sub-tests, type I error rate becomes $1-(1-\alpha)^3 = 0.143$.

k <- 10000

pH0 <- 1:k
pH0.1 <- 1:k
pH0.2 <- 1:k
pH0.3 <- 1:k
n1 <- 4000
n2 <- 4000
n3 <- 4000
n4 <- 4000
n5 <- 4000
n6 <- 4000
p0 <- 0.01

a = 0.05
pD = a

for(i in 1:k){
  x1 <- rbinom(1,n1,p0)
  x2 <- rbinom(1,n2,p0)
  x3 <- rbinom(1,n3,p0)
  x4 <- rbinom(1,n4,p0)
  x5 <- rbinom(1,n5,p0)
  x6 <- rbinom(1,n6,p0)

  pH0.1[i] <- prop.test(c(x1,x2),c(n1,n2), correct = F)$p.value < pD
  pH0.2[i] <- prop.test(c(x3,x4),c(n3,n4), correct = F)$p.value < pD
  pH0.3[i] <- prop.test(c(x5,x6),c(n5,n6), correct = F)$p.value < pD
  pH0[i] <- ifelse((pH0.1[i] | pH0.2[i] | pH0.3[i]) == TRUE, TRUE, FALSE)
}

sum(pH0.1)/k
sum(pH0.2)/k
sum(pH0.3)/k
sum(pH0)/k
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the reply @LuckyPal. The Yate's correction is certainly not necessary in my case, as you say, since I have large samples. My question however, which is obviously very badly phrased, is, how do I do multiple hypothesis testing (which should lead to higher rate of type I error)? What I want to do: multiple hypothesis testing (which should lead to higher rate of type I error unless controlled for). What I am currently doing: Getting 0.05 type I error rate no matter how many (k) tests I do. $\endgroup$ – Pedram Rayat Apr 17 '19 at 7:32
  • $\begingroup$ Your edited version helped @LuckyPal. Thank you for taking the time amigo, I appreciate it a lot since I'm doing my thesis and my supervisor is way too busy to help out. Thanks! $\endgroup$ – Pedram Rayat Apr 18 '19 at 8:08
  • $\begingroup$ that's good to hear @PedramRayat and it's a pity that so often supervisors are too busy to help out. If it indeed answered your question, I would be glad if you could officially accept it as the answer :-) $\endgroup$ – LuckyPal Apr 18 '19 at 8:48
  • 1
    $\begingroup$ It is indeed... Answer accepted! :) @LuckyPal $\endgroup$ – Pedram Rayat Apr 18 '19 at 9:10

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.