# Known alpha and type I error rate simulation

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. • Why do you believe the type I error rate should be "exactly" 5%? – Isabella Ghement Mar 24 '18 at 19:25 • 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. – Peter Flom Mar 24 '18 at 19:58 • 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/…. – whuber Mar 24 '18 at 22:30 • 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. – gung Mar 25 '18 at 0:51 ## 4 Answers 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%

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.

• 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." – whuber Mar 24 '18 at 22:31
• OK, good point. I will amend my answer. – Peter Flom Mar 24 '18 at 22:53

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.

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