library(kableExtra)A Type I error occurs if the null hypothesis is rejected when in fact the null hypothesis is true. An empirical Type I error rate can be computed by a Monte Carlo experiment. The test procedure is replicated a large number of times under the conditions of the null hypothesis. The empirical Type I error rate for the Monte Carlo experiment is the sample proportion of significant test statistics among the replicates. The probability of Type I error is the conditional probability that the null hypothesis is rejected given that H0 is true. Thus, if the test procedure is replicated a large number of times under the conditions of the null hypothesis (H0: \(\ μ=μ0\)), the observed Type I error rate should be at most (approximately) \(\alpha\). The simulation really only investigates empirically whether the method of computing the p-value in t.test (a numerical algorithm) is consistent with the theoretical value \(\alpha\) = 0.05.
The t-test is robust to mild departures from normality. Use Monte Carlo simulation to investigate whether the empirical Type I error rate of the t-test is approximately equal to the nominal significance level \(\alpha\) = 0.05, when the sampled population is non-normal. Discuss the simulation results for the cases where the sampled population is:
set.seed(123456) #keeps the data the same every run
n<- 100 #sample size
alpha <- 0.05 #significance level
mu0 <- 1 #mean of Chi-Squared(1)
m <- 1000 #number of replicates
p <- numeric(m) #storage for p-values
for (j in 1:m) { #Monte Carlo method generating the jth random sample from the null distribution x below
x <- rchisq(n, df = mu0) #calling the distribution under t-test to be of a Chi-Squared Distribution with sample size n and degrees of freedom set to mu0 = 1
ttest <- t.test(x, alternative = "two.sided", mu = mu0) #student two tailed t-test with Chi-Squared(1) distribution with a null hypothesis of mu = mu0
p[j] <- ttest$p.value #computing the proportion of significant tests for each replicate, indexed by j = 1,...,m=1000
}p.hat <- mean(p < alpha) #Type I Error rate
se.hat <- sqrt(p.hat * (1-p.hat) / m) #estimate standard error of Type I Error rate
dt <- round(rbind(p.hat, se.hat),4) #rounding values to 4 decimal places and combining rows of p.hat and se.hat
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate") #labeling the rows
dimnames(dt) <- list(dtnames) #having the data include the names called
dt %>% #creating a neat table using the kableExtra package and giving it a caption
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from Chi-Squared(1) with n = 100") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.0720 |
| Standard Error of Estimate | 0.0082 |
Estimates of Type I error probability will vary, but should be close to the nominal rate \(\alpha\) = 0.05 because all samples were generated under the null hypothesis from the assumed model for a t-test Here, the observed Type I Error rate is 0.0720 and the Standard Error of the estimate is around 0.0082. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.022, when the sampled population is of non-normal \(\chi(1)^2\) and n = 100. Theoretically, the probability of rejecting the null hypothesis of Ho:\(\ μ=μ0\) is 0.0720.
set.seed(123456)
n<- 1000 #increasing sample size to get closer to alpha
alpha <- 0.05
mu0 <- 1
m <- 1000
p <- numeric(m)
for (j in 1:m) {
x <- rchisq(n, df = mu0)
ttest <- t.test(x, alternative = "two.sided", mu = mu0)
p[j] <- ttest$p.value
}p.hat1 <- mean(p < alpha)
se.hat1 <- sqrt(p.hat1 * (1-p.hat1) / m)
dt1 <- round(rbind(p.hat1, se.hat1),4)
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate")
dimnames(dt1) <- list(dtnames)
dt1 %>%
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from Chi-Squared(1) with n = 1,000") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.0640 |
| Standard Error of Estimate | 0.0077 |
Here, the observed Type I Error rate is 0.0640 and the Standard Error of the estimate is around 0.0077. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.014, when the sampled population is of non-normal \(\chi(1)^2\) and n = 1,000. Theoretically, the probability of rejecting the null hypothesis of Ho: \(\ μ=μ0\) is 0.0640 which is 0.008 closer than with n = 100.
set.seed(123456)
n<- 10000 #increasing sample size to get closer to alpha
alpha <- 0.05
mu0 <- 1
m <- 1000
p <- numeric(m)
for (j in 1:m) {
x <- rchisq(n, df = mu0)
ttest <- t.test(x, alternative = "two.sided", mu = mu0)
p[j] <- ttest$p.value
}p.hat2 <- mean(p < alpha)
se.hat2 <- sqrt(p.hat2 * (1-p.hat2) / m)
dt2 <- round(rbind(p.hat2, se.hat2),4)
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate")
dimnames(dt2) <- list(dtnames)
dt2 %>%
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from Chi-Squared(1) with n = 10,000") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.0460 |
| Standard Error of Estimate | 0.0066 |
Here, the observed Type I Error rate is 0.0460 and the Standard Error of the estimate is around 0.0066. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.004, when the sampled population is of non-normal \(\chi(1)^2\) and n = 10,000. Theoretically, the probability of rejecting the null hypothesis of Ho: \(\ μ=μ0\) is 0.0460 which is 0.010 closer than with n = 1,000.
set.seed(123456) #keeps the data the same every run
n<- 100 #sample size
alpha <- 0.05 #significance level
mu0 <- 1 #mean of U(0,2)
m <- 1000 #number of replicates
p <- numeric(m) #storage for p-values
for (j in 1:m) { #Monte Carlo method generating the jth random sample from the null distribution x below
x <- runif(n,0,2) #calling the distribution under t-test to be of a Uniform(0,2) Distribution with sample size n and degrees of freedom set to mu0 = 1
ttest <- t.test(x, alternative = "two.sided", mu = mu0) #student two tailed t-test with U(0,2) distribution with a null hypothesis of mu = mu0
p[j] <- ttest$p.value #computing the proportion of significant tests for each replicate, indexed by j = 1,...,m=1000
}p.hat3 <- mean(p < alpha) #Type I Error rate
se.hat3 <- sqrt(p.hat3 * (1-p.hat3) / m) #estimate standard error of Type I Error rate
dt3 <- round(rbind(p.hat3, se.hat3),4) #rounding values to 4 decimal places and combining rows of p.hat and se.hat
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate") #labeling the rows
dimnames(dt3) <- list(dtnames) #having the data include the names called
dt3 %>% #creating a neat table using the kableExtra package and giving it a caption
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from U(0,2) with n = 100") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.0440 |
| Standard Error of Estimate | 0.0065 |
Here, the observed Type I Error rate is 0.0440 and the Standard Error of the estimate is around 0.0065. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.006, when the sampled population is of non-normal \(\ U(0,2)\) dist. and n = 100. Theoretically, the probability of rejecting the null hypothesis of Ho: \(\ μ=μ0\) is 0.0440.
set.seed(123456)
n<- 1000 #increasing sample size to get closer to alpha
alpha <- 0.05
mu0 <- 1
m <- 1000
p <- numeric(m)
for (j in 1:m) {
x <- runif(n,0,2)
ttest <- t.test(x, alternative = "two.sided", mu = mu0)
p[j] <- ttest$p.value
}p.hat4 <- mean(p < alpha)
se.hat4 <- sqrt(p.hat4 * (1-p.hat4) / m)
dt4 <- round(rbind(p.hat4, se.hat4),4)
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate")
dimnames(dt4) <- list(dtnames)
dt4 %>%
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from U(0,2) with n = 1,000") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.038 |
| Standard Error of Estimate | 0.006 |
Here, the observed Type I Error rate is 0.0380 and the Standard Error of the estimate is around 0.0006. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.012, when the sampled population is of non-normal \(\ U(0,2)\) dist. and n = 1,000. Theoretically, the probability of rejecting the null hypothesis of Ho: \(\ μ=μ0\) is 0.0380 which is actually 0.006 further away than with n = 100.
set.seed(123456)
n<- 10000 #increasing sample size to get closer to alpha
alpha <- 0.05
mu0 <- 1
m <- 1000
p <- numeric(m)
for (j in 1:m) {
x <- runif(n,0,2)
ttest <- t.test(x, alternative = "two.sided", mu = mu0)
p[j] <- ttest$p.value
}p.hat5 <- mean(p < alpha)
se.hat5 <- sqrt(p.hat5 * (1-p.hat5) / m)
dt5 <- round(rbind(p.hat5, se.hat5),4)
dtnames <- c("Empirical Type-I Error", "Standard Error of Estimate")
dimnames(dt5) <- list(dtnames)
dt5 %>%
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Population from U(0,2) with n = 10,000") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | 0.051 |
| Standard Error of Estimate | 0.007 |
Here, the observed Type I Error rate is 0.0510 and the Standard Error of the estimate is around 0.007. As you can see, the Empirical Type I Error rate differs from the nominal significance level \(\alpha\) = 0.05, by approximately 0.001, when the sampled population is of non-normal \(\ U(0,2)\) dist. and n = 100. Theoretically, the probability of rejecting the null hypothesis of Ho: \(\ μ=μ0\) is 0.0510 which is 0.013 closer than with n = 1,000.
a <- c(p.hat, p.hat1, p.hat2, p.hat3, p.hat4, p.hat5, se.hat, se.hat1, se.hat2, se.hat3, se.hat4, se.hat5) #creating a vector of all of the different Type I Error rates and standard errors
tab <- matrix(rep(a, times = 1), ncol=2, nrow = 6) #creating a matrix of the values
colnames(tab) <- c("Empirical Type-I Error", "Standard Error of Estimate") #labeling the columns
rownames(tab) <- c("Chi-Squared(1): n = 100", "n = 1,000", "n = 10,000", "U(0,2): n = 100", "n = 1,000", "n = 10,000") #labeling the rows
tab <- as.table(round(tab, digits = 4)) #creating the data table of the values combined and rounding to 4 decimal places
tab %>%
kbl(caption = "Monte Carlo Simulation with Non-Normal Sampled Populations") %>%
kable_classic(full_width = F, html_font = "Cambria")| Empirical Type-I Error | Standard Error of Estimate | |
|---|---|---|
| Chi-Squared(1): n = 100 | 0.072 | 0.0082 |
| n = 1,000 | 0.064 | 0.0077 |
| n = 10,000 | 0.046 | 0.0066 |
| U(0,2): n = 100 | 0.044 | 0.0065 |
| n = 1,000 | 0.038 | 0.0060 |
| n = 10,000 | 0.051 | 0.0070 |
As you can see from the table, the t-test is robust to mild departures from normality. Thus, when the test procedure was replicated a large number of times under the conditions of the null hypothesis Ho: \(\ μ=μ0\), the observed Type I error rate got closer to \(\alpha\) = 0.05 as the sample size increased in both non-normal distributions.