Homework Assignment

Design a small simulation study to investigate the empirical performance of the two tests mentioned during exercise class. Firstly, simulate the data from the null hypothesis and perform both tests. Find out, what is the empirical probability of the first type error (i.e., the overall proportion of rejecting the null hypothesis, if the null hypothesis holds). Consider some data from the alternative hypothesis and investigate what is the empirical power of the test (i.e., the empirical probability of not rejecting the null, if the null hypothesis does not hold). Use different sample size and different alternative hypothesis to get some meaningful insight into the performance of both tests. Provide some graphical outputs and post the figures (together with some short comments).

At first we simulate data from the bivariate normal distribution \(N_{2}\Big(\Big(\begin{array}{c}1\\1\end{array}\Big), \Big( \begin{array}{cc}1 & 0.5\\0.5 & 2 \end{array} \Big)\Big)\), and test the null hypothesis \[ H_{0}: \mu_{1} - \mu_2 = 0.\] We consider two situations: either the variance-covariance matrix \(\Sigma\) is known, or \(\Sigma\) is unknown. We use the test statistics derived during exercise class for both cases. We choose the sample size equal to 10,100,1000 and 10000 and compute 500 times the corresponding tests. We compute the proportion of rejected tests and plot the empirical probability of the first type of error.

library("mvtnorm")
set.seed(1234)
samplesz <- c(10,100,1000,10000)
mu <-  c(1, 1)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
rej_S_known <- c(rep(0,length(samplesz)))
rej_S_unknown <- c(rep(0,length(samplesz)))
a <- 0
A <- c(1, -1) 
#n_tests <- 1000
n_tests <- 500
for (n in samplesz){
for(i in 1:n_tests){
sample <- rmvnorm(n, mu, Sigma)
###compute test for Sigma known
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
#if we reject that add 1 to rej_S_unknown
if (Tvalue>qchisq(0.95,1)) rej_S_known[log10(n)]=rej_S_known[log10(n)]+1 
#compute test for Sigma unknown
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
if (TTvalue>qf(0.95,1,99)) rej_S_unknown[log10(n)]=rej_S_unknown[log10(n)]+1 
}}
sig_S_known <-rej_S_known/n_tests
sig_S_unknown <-rej_S_unknown/n_tests
plot(c(1:4),sig_S_known, pch=21, bg="deepskyblue3", ylim = c(0,0.2), xlab = "Log_10(n)", ylab = "The empirical probability of the first type error", main = "The empirical probability of the first type of error")
lines(c(1:4),sig_S_known, col = "deepskyblue1", lwd = 2)   
points(c(1:4),sig_S_unknown, pch=21, bg="firebrick4")
lines(c(1:4),sig_S_unknown, col="firebrick1", lwd = 2)
abline(h=0.05, col="black", lty=2)
legend(3, 0.2, legend = c("Sigma known","Sigma unknown"), title = "Type of test", col = c("deepskyblue1","firebrick1"), lty = 1, lwd = 2)

From the plot we observe that both test (for (un)known sigma matrix) have the expected probability of the first type of error around 0.05 and the lines do not differ to much for the two types of test.

Now we simulate the data from the bivariate normal distribution \(N_{2}\Big(\Big(\begin{array}{c}1\\1.1\end{array}\Big), \Big( \begin{array}{cc}1 & 0.5\\0.5 & 2 \end{array} \Big)\Big)\), and test the null hypothesis \[ H_{0}: \mu_{1} - \mu_2 = 0.\]. So we simulate the data for alternative and look for power of the test again for both two versions for (un)known sigma matrix.

library("mvtnorm")
set.seed(1234)
samplesz <- c(10,100,1000,10000,100000)
mu <-  c(1, 1.1)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
rej_S_known <- c(rep(0,length(samplesz)))
rej_S_unknown <- c(rep(0,length(samplesz)))
a <- 0
A <- c(1, -1) 

n_tests <- 100
for (n in samplesz){
for(i in 1:n_tests){
sample <- rmvnorm(n, mu, Sigma)
###compute test for Sigma known
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
#if we reject that add 1 to rej_S_unknown
if (Tvalue>qchisq(0.95,1)) rej_S_known[log10(n)]=rej_S_known[log10(n)]+1 
#compute test for Sigma unknown
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
if (TTvalue>qf(0.95,1,99)) rej_S_unknown[log10(n)]=rej_S_unknown[log10(n)]+1 
}}
sig_S_known <-rej_S_known/n_tests
sig_S_unknown <-rej_S_unknown/n_tests
plot(c(1:length(samplesz)),sig_S_known, pch=21, bg="deepskyblue3", ylim = c(0,1.2), xlab = "Log_10(n)", ylab = "The empirical power of the test", main = "Power of test, data for mu changed by 0.1")
lines(c(1:length(samplesz)),sig_S_known, col = "deepskyblue1", lwd = 2)   
points(c(1:length(samplesz)),sig_S_unknown, pch=21, bg="firebrick4")
lines(c(1:length(samplesz)),sig_S_unknown, col="firebrick1", lwd = 2)
abline(h=1, col="black", lty=2)
legend(1, 1.2, legend = c("Sigma known","Sigma unknown"), title = "Type of test", col = c("deepskyblue1","firebrick1"), lty = 1, lwd = 2)

From the plots we observe that with increasing sample size the power of the test goes to one which is fine. The curves for both type of test do not differ too much.

We also try to simulate data under different violation of null hypothesis. We simulate the data from multivariate t-distribution with two degrees of freedom mean vector \((1,1)^\top\) and sigma matrix \(\Big( \begin{array}{cc}1 & 0.5\\0.5 & 2 \end{array} \Big)\) (for details see function draw.d.variate.t) so we only violated the assumption of multivariate normal distribution. We plot the again the power of test for both versions.

library("mvtnorm")
library("MultiRNG")
set.seed(1234)
samplesz <- c(10,100,1000,10000,100000)
mu <-  c(1, 1)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
rej_S_known <- c(rep(0,length(samplesz)))
rej_S_unknown <- c(rep(0,length(samplesz)))
a <- 0
A <- c(1, -1) 
#n_tests <- 1000
n_tests <- 100
for (n in samplesz){
for(i in 1:n_tests){
sample <- draw.d.variate.t(2,n,2,mu,Sigma)
###compute test for Sigma known
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
#if we reject that add 1 to rej_S_unknown
if (Tvalue>qchisq(0.95,1)) rej_S_known[log10(n)]=rej_S_known[log10(n)]+1 
#compute test for Sigma unknown
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
if (TTvalue>qf(0.95,1,99)) rej_S_unknown[log10(n)]=rej_S_unknown[log10(n)]+1 
}}
sig_S_known <-rej_S_known/n_tests
sig_S_unknown <-rej_S_unknown/n_tests
plot(c(1:length(samplesz)),sig_S_known, pch=21, bg="deepskyblue3", ylim = c(0,1.2), xlab = "Log_10(n)", ylab = "The empirical power of the test", main = "Power of test, data from t-distribution")
lines(c(1:length(samplesz)),sig_S_known, col = "deepskyblue1", lwd = 2)   
points(c(1:length(samplesz)),sig_S_unknown, pch=21, bg="firebrick4")
lines(c(1:length(samplesz)),sig_S_unknown, col="firebrick1", lwd = 2)
abline(h=1, col="black", lty=2)
abline(h=0.05, col="black", lty=2)
legend(1, 1.2, legend = c("Sigma known","Sigma unknown"), title = "Type of test", col = c("deepskyblue1","firebrick1"), lty = 1, lwd = 2)

We observe that the power of the test is low and do not improve with the increase of the sample size. We can see small difference between known and unknown sigma matrix (with higher power of the test for known sigma matrix).