The data

First of all, let’s simulate some data and see if the means and the variance correspond. The data follow a multivariate normal distribution such that : \[{X \sim\mathcal{N}_2(\begin{pmatrix}1\\2\end{pmatrix},\begin{pmatrix}1 & 0.5\\0.5 & 1\end{pmatrix})}\].

n=200
mu <- c(1,2)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
set.seed(1234)
sample <- rmvnorm(n, mu, Sigma)

sample=as.data.frame(sample)
sapply(sample, mean, na.rm = TRUE)                  ### not parallel
##       V1       V2 
## 1.012694 2.007053
cov(sample)
##           V1        V2
## V1 1.0396159 0.6051607
## V2 0.6051607 2.1386796

The tests

Now we will run two test about the hypothesis : \[ \mathcal{H}_0 : \mathcal{A}\mu = a \] vs \[ \mathcal{H}_1 : \mathcal{A}\mu \ne a \] where \(\mathcal{A} =\begin{pmatrix}2\\-1\end{pmatrix}\) and \(a = 0\).

With \(\Sigma\) known

We know that, under \(\mathcal{H}_0\), the likelihood ratio statistic follows a chi2 statistic, with one degree of freedom.

We calculate it, such as the quantile of the chi2 distribution, with a level of 95 %.

a <- 0
A <- c(2, -1) 
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
Tvalue
##            [,1]
## [1,] 0.01680831
qchisq(0.95,df=1)
## [1] 3.841459

Because the T value is way under the quantile, we can not reject \(\mathcal{H}_0\), therefore we accept it.

With \(\Sigma\) unknown

Although the process is quite similar, we now perform a test where the statistic follows a Fisher distribution with 1 and (n-1) degrees of freedom.

TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
TTvalue
##            [,1]
## [1,] 0.01725708
qf(0.95,df1 = 1,df2=(n-1))
## [1] 3.888613

Once again, we can not reject \(\mathcal{H}_0\), because the T value is way under the quantile of the Fisher distribution.

Type errors

We can find the type erros of these tests :

#install.packages('pwr')
library(pwr)
## Warning: le package 'pwr' a été compilé avec la version R 4.1.3
#?pwr.chisq.test
pwr.chisq.test(w=Tvalue,N=n,df=1,sig.level = 0.95)
## 
##      Chi squared power calculation 
## 
##               w = 0.01680831
##               N = 200
##              df = 1
##       sig.level = 0.95
##           power = 0.951391
## 
## NOTE: N is the number of observations
#?pwr.f2.test
pwr.f2.test(u=1,v=(n-1),f2=TTvalue,sig.level=0.95)
## 
##      Multiple regression power calculation 
## 
##               u = 1
##               v = 199
##              f2 = 0.01725708
##       sig.level = 0.95
##           power = 0.9911539

We know having a high power number is a good criteria to decide if the test is appropriate or not. Thus, both tests look interesting.

We can calculate the Type Error II, which is : \[ Type Error II = 1 - power\].

Therefore, the empirical probability of not rejecting the null, if the null hypothesis does not hold is in each test: \(1 - 0.95 = 0.05\) .

Data from the alternative hypothesis

Now, we will look into the same tests, but with data that does nos follow the null hypothesis such that : \[ \mathcal{H}_1 : \mathcal{A}\mu \ne a \]

We will run the same two tests with the hypothesis of \(\Sigma\) know first and then unknown.

#new data
n=200
mu <- c(1,2)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
set.seed(1234)
sample <- rmvnorm(n, mu, Sigma)
sample=as.data.frame(sample)
#sapply(sample, mean, na.rm = TRUE)                  

#Sigma known
a <- 2
A <- c(2,-1) 
Tvalue <- n * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
Tvalue
##          [,1]
## [1,] 196.3498
qchisq(0.95,df=1)
## [1] 3.841459
#Sigma unknown
TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
TTvalue
##          [,1]
## [1,] 201.5922
qf(0.95,df1 = 1,df2=(n-1))
## [1] 3.888613

Here, we reject the null hypothesis in both cases because the statistic is way bigger than the distribution quantile.

Global Perfomance

Let’s iterate the test using n to see how the test behave while changing the n.

mu <- c(1,2)
Sigma <- matrix(c(1, 0.5, 0.5, 2),2,2)
a <- 0
A <- c(2, -1) 

first_test = NULL
second_test = NULL
type_error_2_1 = NULL
type_error_2_2 = NULL

for (i in 10:5000){
                 n = i 
                 sample = rmvnorm(n, mean = mu, sigma = Sigma)
                 
                 #First test
                 Tvalue = n* t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% Sigma %*% A) %*% (A %*% apply(sample, 2, mean) - a)
                 first_test = c(first_test, Tvalue > qchisq(0.95, df = 1))
                 type_error_2_1 = c(type_error_2_1,1 - pwr.chisq.test(w=Tvalue,N=n,df=1,sig.level = 0.95)$power)
                 
                 #Second test
                 TTvalue <- (n - 1) * t(A %*% apply(sample, 2, mean) - a) %*% solve(A %*% cov(sample) %*% A) %*% (A %*% apply(sample, 2, mean) - a)
                 second_test = c(second_test, TTvalue > qf(0.95, df1 = 1, df2 = n-1))
                 type_error_2_2 = c(type_error_2_2, 1 - pwr.f2.test(u = 1,v = (n-1), f2 = TTvalue, sig.level = 0.05)$power)
}

mean(type_error_2_1)
## [1] 0.00662695
mean(type_error_2_2)
## [1] 0.04201238
sum(first_test)/4990
## [1] 0.05490982
sum(second_test)/4990
## [1] 0.05410822

As we can see, while raising \(n\), our tests remain performant. We have means about 0.05 on 4990 tests and when we sum up all the cases where we reject the null hypothesis and divide them by the number of test, we obtain the small probability of 0.05.

Graphical outputs

if we plot the Type Error II for both tests, we obtain these graphics.

boxplot(type_error_2_1)

boxplot(type_error_2_2)

We can see most of the values are close to zero, which is a good thing. However, while the first boxplot has only values between \(0\) and \(0.05\), the second has values rising to more than \(0.8\).

This implies that if the \(\Sigma\) is known, the tests are more precise.