Known alpha and type I error rate simulation

Question

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.

Answer 1

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%

Answer 2

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:

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.

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

Answer 3

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.

Answer 4

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.