2
$\begingroup$

I am not sure why running the following code in R does not result in a type I error rate of exactly 5% given alpha = 0.05 when set.seed is set to 1. Is there a reason why running the simulation below gives a type I error rate of 4% when it should be exactly 5% using set.seed(1)? Sometimes when you run the simulation, you get a p-value <0.05 in 3.2% of the simulations, is this due to rnorm drawing a random normal number each time?

set.seed(1)
n=1000 # testing 1,000 times
t1err=0
for (i in 1:n){
  x=rnorm(100, 0, 1)
  if (((t.test(x, mu=0))$p.value)<=0.05) (t1err=t1err+1) 
}
cat("Type I error rate in percentage is", (t1err/n)*100,"%")

Why is it that seed = 123, p = 3.2% sometimes with n = 1000? This seems far off from the stated true value of 5%? Is this due to how R draws random normal numbers each time using rnorm regardless of set.seed? Supposedly the function of set.seed is to make the simulations more reproducible each time.

$\endgroup$
  • 4
    $\begingroup$ Why do you believe the type I error rate should be "exactly" 5%? $\endgroup$ – Isabella Ghement Mar 24 '18 at 19:25
  • 2
    $\begingroup$ This is asked about R, but the same thing would happen with other simulators. It's not really about R, so I vote to leave it open. $\endgroup$ – Peter Flom Mar 24 '18 at 19:58
  • $\begingroup$ When you are patient enough, you can create non-random sequences with any pseudorandom number generator. The technique is explained and illustrated at stats.stackexchange.com/questions/38063/…. $\endgroup$ – whuber Mar 24 '18 at 22:30
  • $\begingroup$ Please register &/or merge your accounts (you can find information on how to do this in the My Account section of our help center), then you will be able to edit & comment on your own question. $\endgroup$ – gung Mar 25 '18 at 0:51
2
$\begingroup$

If you simulate a lot more tests it gets closer to 5% - I think this is a case of a slower convergence to the "true" value than you seem to expect:

    set.seed(1)
    n=100000 # testing 100,000 times
    t1err=0
    for (i in 1:n){
      x=rnorm(100, 0, 1)
      if (((t.test(x, mu=0))$p.value)<=0.05) (t1err=t1err+1) 
    }
    cat("Type I error rate in percentage is", (t1err/n)*100,"%")

gives 4.894%

$\endgroup$
1
$\begingroup$

It's because rnorm is random draws. If you use a different seed, you will get different values.

seed = 12, p = 4.9%
seed = 123, p = 3.2%
seed = 1234, p = 5.4%

and, as @dlmoore pointed out, 1,000 is not very many. If you use 100,000, then with

seed = 1, p = 4.89%
seed = 12, p = 4.89%
seed = 123, p = 5.04%
seed = 1234, p = 5.13%

but there is no reason for it to equal exactly 5.00000.

$\endgroup$
  • 2
    $\begingroup$ Be careful with "never," Peter: in this Binomial simulation of 1000 draws there is almost a 5.8% chance of obtaining exactly a 5% result. In a simulation 100 times longer (100,000) the chance drops to a tenth of that, but a value greater then a half percent is still appreciably different from "never." $\endgroup$ – whuber Mar 24 '18 at 22:31
  • 1
    $\begingroup$ OK, good point. I will amend my answer. $\endgroup$ – Peter Flom Mar 24 '18 at 22:53
1
$\begingroup$

Remember that p-values, test decisions, test statistics etc. are (pseudo-)random variables. You have some realization off one for each (simulated) dataset, but the expected value is solely what you expect in the long-run.

In fact, here it is clear that for a single simulated dataset the estimated type 1 error rate would soothe be 0 or 1. For two datasets 0, 0.5 or 1 and so on. You first have any possibility of exactly 0.05 for 20 simulated datasets.

If you consider the new random variable that is the estimated type 1 error rate across a finite number $m$ of datasets, its expectation is 0.05 for all values of $m>0$ in your example, but it's variance decreases with $m$, but even with rounding to 4 decimals it takes a while until your are almost certain to get 0.0500.

$\endgroup$
0
$\begingroup$

Basically you are flipping a coin 1000 times and it will fall on heads (error) 5% of the time. The distribution for such occurrences is a binomial distribution:

distribution for number of errors

So every time you redo the 1000 experiments you might observe a different number of failures.

  • The expected number of failures will be 5% (or 50 occurrences in 1000), when the null hypothesis is true.
  • but the observed number of failures is a variable which will be distributed according to the above image and is not necessarily equal to 5% every time.

The p-value is, in this case, just a (convenient) re-parameterization of the t-value, which is a random variable


See also below an example of your simulation repeated for 10 000 different seeds. You will get a different number of errors each time. However, the distribution of these different numbers resembles the distribution in the image above.

example of distribution for different seeds

errors <- sapply(1:10000, FUN = function(x) {
  set.seed(x)
  n=1000 # testing 1,000 times
  t1err=0
  for (i in 1:n){
    x=rnorm(100, 0, 1)
    if (((t.test(x, mu=0))$p.value)<=0.05) (t1err=t1err+1) 
  }
  t1err
})

layout(matrix(c(1,2), 1, 2, byrow = TRUE),width=c(4,1))
plot(1:1000,errors[1:1000], 
     xlab="used seed", ylab="number of errors",
     ylim=c(25,75),
     pch=21,cex=0.7,
     col=adjustcolor("black",alpha.f=1),bg=adjustcolor("black",alpha.f=1))
title("number of errors for different seeds  \n (1 000 shown)")

yhist <- hist(errors, plot = FALSE,breaks=seq(25,75,by=2.5))
margins <- par()$mar
par(mar = c(margins[1],0,margins[3],1))
barplot(yhist$counts, axes = 0, horiz = TRUE)
title("histogram  \n (for 10 000)")
$\endgroup$

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.