two_var_matrices_test_G2 <- function(s, n) {
m <- nrow(s[[1]]) # number of variables
N1 <- n[[1]] # sample size 1
N2 <- n[[2]] # sample size 2
n1 <- N1 - 1 # Using expression (1.2)
n2 <- N2 - 1 # Using expression (1.2)
n <- n1 + n2
# Matrices V1, V2 and V1_plus
V1 <- s[[1]] * n1 # Using expression (1.2)
V2 <- s[[2]] * n2 # Using expression (1.2)
V1_plus <- pracma::pinv(V1)
# Obtaining V
V <- V1 + V2
# My trace function
my_trace <- function(x) sum(diag(x))
# Estimating a1, a2 and b
a1_hat <- my_trace(V) / (n*m)
denomi <- (n-1)*(n+2)*m
a2_hat <- (my_trace(V %*% V) + my_trace(V)^2 / n) / denomi
b_hat <- a1_hat^2 / a2_hat
# The statistic and the p-value
G2 <- m * b_hat * my_trace(V1_plus %*% V2)
p.value <- pchisq(q=G2, df=n1*n2, lower.tail=FALSE)
# The usual way to report a test in R is:
method <- "G2 test for homogeneity of covariances"
statistic <- G2
names(statistic) <- c("G2")
parameter <- c(n1*n2)
names(parameter) <- c("df")
alternative <- "the covariance matrices are different \n"
estimate <- "Due to the high value of m, matrices S1 and S2 are not displayed."
data.name <- "this test uses summarized data"
result <- list(statistic = statistic,
parameter = parameter,
p.value = p.value,
estimate = estimate,
alternative = alternative,
method = method,
data.name = data.name)
class(result) <- "htest"
return(result)
}1 Introducción
En esta publicación vamos a implementar en R las pruebas propuestas en el artículo de Srivastava y Yanagihara (2010).
2 Utilidad de estas pruebas
Las pruebas mostradas en esta aplicación son útiles cuando el número de variables \(m\) excede el número de observaciones \(n_i\) tomada de cada una de las dos poblaciones. Un ejemplo donde sucede con frecuencia que se tienen más variables que observaciones es en estudios genéticos.
El objetivo de las pruebas se puede resumir en el siguiente conjunto de hipótesis:
\[\begin{equation} \label{eq1} \begin{split} H_0 &: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2 \\ H_1 &: \boldsymbol{\Sigma}_1 \neq \boldsymbol{\Sigma}_2 \end{split} \end{equation}\]
3 Prueba \(G_2\)
El estadístico de esta prueba está dado por
\[ G_2 = m \, \hat{b} \, \text{tr}(V_1^+ V_2), \]
donde \(m\) corresponde al número de variables. Las restantes cantidades están definidas en el artículo Srivastava y Yanagihara (2010). La distribución del estadístico \(G_2\) bajo el supuesto de \(H_0\) verdadera es \(\chi^2_{n_1 \times n_2}\) siempre que \(m \rightarrow \infty\).
A continuación se presenta la función de R two_var_matrices_test_G2 creada para realizar la prueba. A la función se le debe ingresar una lista con las matrices de varianzas y covarianzas por medio del argumento s y una lista con los tamaños muestrales por medio del argumento n.
3.1 Ejemplo 1 y \(H_0\) verdadera
En este primer ejemplo vamos a simular dos muestras de dos poblaciones normales \(N_m(\boldsymbol{0}, \boldsymbol{\Sigma}_i)\) con \(m=7\). La primer muestra tendrá \(n_1=4\) observaciones y la segunda \(n_2=5\) observaciones. En este ejemplo el valor de \(m\) no es muy grande, esto se hace solo como ilustración.
Ambas poblaciones tienen vector de medias \(\boldsymbol{0}\). La matriz \(\boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\) serán matrices diagonal con todas las varianzas iguales a 1. Esto garantiza que la hipótesis nula \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\) sea verdadera y lo que esperamos es que la prueba sea capaz de detectar esta situación y que no rechace \(H_0\).
En el código de abajo se definen los valores, se fijan las semillas para que el lector pueda obtener los mismos resultados y se crean las listas s y n.
library(mvtnorm)
m <- 7
n1 <- 4
n2 <- 5
set.seed(123)
x1 <- rmvnorm(n=n1, mean=rep(0, m), sigma=diag(m))
s1 <- var(x1)
set.seed(123)
x2 <- rmvnorm(n=n2, mean=rep(0, m), sigma=diag(m))
s2 <- var(x2)
s <- list(s1, s2)
n <- list(n1, n2)Ahora vamos a utilizar la función creada de la siguiente manera.
res1 <- two_var_matrices_test_G2(s=s, n=n)
res1
G2 test for homogeneity of covariances
data: this test uses summarized data
G2 = 6.6802, df = 12, p-value = 0.878
alternative hypothesis: the covariance matrices are different
sample estimates:
[1] "Due to the high value of m, matrices S1 and S2 are not displayed."De la salida anterior vemos que el valor-P es mayor que el nivel de significancia usual del 5%, esto significa que no hay evidencias para rechazar \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\), como era de esperarse.
A continuación vamos a ilustrar el valor-P encontrado usando la distribución del estadístico.
curve(dchisq(x, df=res1$parameter), from=0, to=40,
lwd=3, xlab="Statistic")
segments(x0=res1$statistic, y0=0,
x1=res1$statistic, y1=dchisq(x=res1$statistic, df=res1$parameter),
lty="dashed", col="tomato", lwd=2)
Como el estadístico es 6.6801706, en la figura se observa una línea de color naranja en este valor, el área a la derecha de este valor y debajo de la curva es 0.8780056.
3.2 Ejemplo 2 y \(H_0\) falsa
En este primer ejemplo vamos a simular dos muestras de dos poblaciones normales \(N_m(\boldsymbol{0}, \boldsymbol{\Sigma}_i)\) con \(m=7\). La primer muestra tendrá \(n_1=4\) observaciones y la segunda \(n_2=5\) observaciones. En este ejemplo el valor de \(m\) no es muy grande, esto se hace solo como ilustración.
Ambas poblaciones tienen vector de medias \(\boldsymbol{0}\). La matriz \(\boldsymbol{\Sigma}_1\) será una matriz diagonal con todas las varianzas iguales a 1, mientras que la matriz \(\boldsymbol{\Sigma}_2\) será una matriz diagonal con varianzas iguales a 4. Esto garantiza que la hipótesis nula \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\) sea verdadera y lo que esperamos es que la prueba sea capaz de detectar esta situación y que rechace \(H_0\).
En el código de abajo se definen los valores, se fijan las semillas para que el lector pueda obtener los mismos resultados y se crean las listas s y n.
library(mvtnorm)
m <- 7
n1 <- 4
n2 <- 5
set.seed(123)
x1 <- rmvnorm(n=n1, mean=rep(0, m), sigma=diag(m))
s1 <- var(x1)
set.seed(123)
x2 <- rmvnorm(n=n2, mean=rep(0, m), sigma=4*diag(m))
s2 <- var(x2)
s <- list(s1, s2)
n <- list(n1, n2)Ahora vamos a utilizar la función creada de la siguiente manera.
res2 <- two_var_matrices_test_G2(s=s, n=n)
res2
G2 test for homogeneity of covariances
data: this test uses summarized data
G2 = 27.309, df = 12, p-value = 0.006974
alternative hypothesis: the covariance matrices are different
sample estimates:
[1] "Due to the high value of m, matrices S1 and S2 are not displayed."De la salida anterior vemos que el valor-P es menor que el nivel de significancia usual del 5%, esto significa que hay evidencias para rechazar \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\), como era de esperarse.
A continuación vamos a ilustrar el valor-P encontrado usando la distribución del estadístico.
curve(dchisq(x, df=res2$parameter), from=0, to=40,
lwd=3, xlab="Statistic")
segments(x0=res2$statistic, y0=0,
x1=res2$statistic, y1=dchisq(x=res2$statistic, df=res2$parameter),
lty="dashed", col="tomato", lwd=2)
Como el estadístico es 27.3087735, en la figura se observa una línea de color naranja en este valor, el área a la derecha de este valor y debajo de la curva es 0.0069736.
3.3 Ejemplo 3 con \(m >>> n_i\)
En este ejemplo vamos a simular dos muestras de dos poblaciones normales (\(N_m(\boldsymbol{0}, \boldsymbol{\Sigma}_i)\)) con \(m=30\). La primer muestra tendrá \(n_1=4\) observaciones y la segunda \(n_2=5\) observaciones. En este ejemplo el valor de \(m=30\) cumple con la condición principal de las pruebas, mayor dimensión que número de observaciones.
Ambas poblaciones tienen vector de medias \(\boldsymbol{0}\). La matriz \(\boldsymbol{\Sigma}_1\) será una matriz diagonal con todas las varianzas iguales a 1, mientras que la matriz \(\boldsymbol{\Sigma}_2\) será una matriz diagonal con varianzas iguales a 15. Esto garantiza que la hipótesis nula \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\) no sea verdadera y lo que esperamos es que la prueba sea capaz de detectar esta situación.
En el código de abajo se definen los valores, se fijan las semillas para que el lector pueda obtener los mismos resultados y se crean las listas s y n.
library(mvtnorm)
m <- 30
n1 <- 4
n2 <- 5
set.seed(123)
x1 <- rmvnorm(n=n1, mean=rep(0, m), sigma=diag(m))
s1 <- var(x1)
set.seed(123)
x2 <- rmvnorm(n=n2, mean=rep(0, m), sigma=15*diag(m))
s2 <- var(x2)
s <- list(s1, s2)
n <- list(n1, n2)Ahora vamos a utilizar la función creada de la siguiente manera.
res3 <- two_var_matrices_test_G2(s=s, n=n)
res3
G2 test for homogeneity of covariances
data: this test uses summarized data
G2 = 122.9, df = 12, p-value < 2.2e-16
alternative hypothesis: the covariance matrices are different
sample estimates:
[1] "Due to the high value of m, matrices S1 and S2 are not displayed."De la salida anterior vemos que el valor-P es menor que el nivel de significancia usual del 5%, esto significa que hay evidencias para rechazar \(H_0: \boldsymbol{\Sigma}_1 = \boldsymbol{\Sigma}_2\), como era de esperarse.
A continuación vamos a ilustrar el valor-P encontrado usando la distribución del estadístico.
curve(dchisq(x, df=res3$parameter), from=0, to=60,
lwd=3, xlab="Statistic")
segments(x0=res3$statistic, y0=0,
x1=res3$statistic, y1=dchisq(x=res3$statistic, df=res3$parameter),
lty="dashed", col="tomato", lwd=2)
Como el estadístico es 122.8996297, en la figura se observa una línea de color naranja en este valor, el área a la derecha de este valor y debajo de la curva es 1.6304521^{-20}.
4 Prueba \(J_2\)
El estadístico de esta prueba está dado por
\[ J_2 = bla bla bla, \]
5 Prueba \(T_2\)
El estadístico de esta prueba está dado por
\[ T_2 = bla bla bla, \]
6 Prueba \(Q_2\)
El estadístico de esta prueba está dado por
\[ Q_2 = bla bla bla, \]