Introducción

Para implementar la prueba de Johansen de igual de los vectores de medias de dos poblaciones normales es posible urilizar el paquete stests de \(R\). En este paquete se pueden encontrar otras pruebas útiles en el análisis multivariado. Para instalar y utilizar este paquete se escriben las siguentes lineas de código:

if (!require('devtools')) install.packages('devtools')
devtools::install_github('fhernanb/stests', force=TRUE)
library(stests)

Prueba de Johansen

Sea \(\bf{X_{i1}},\,\bf{X_{i2}},\,\ldots,\,\bf{X_{in_i}},\,\) una muestra aleatoria con una distribución normal \(p\)-variada \(N(\bf\mu_i,\Sigma_i)\) e \(i=1,\,2\). Para estas dos poblaciones se desean contrastar las siguientes dos hipótesis: \[H_0:\bf\mu_1=\bf\mu_2\,\, vs\, \, \bf\mu_1\neq\bf\mu_1 \] Este contraste es conocido como la el problema de Behrens-Fisher. Para resolver este problema se pueden utilizar diferentes Test estadísticos. En esta publicación se centra la atención en la Prueba de Johansen (1980). En este caso el estadístico de prueba está dado por:

\[T^2= \left(\bf\bar{X}_1-\bf\bar{X}_2\right)^\top \bf{\tilde{S}}^{-1} \left(\bf\bar{X}_1-\bf\bar{X}_2\right)\] donde \(\bf\bar{X}_i\) es la media muestral de la \(i-\)ésima población para \(i=1,\,2\), y \(\tilde{S}\) está dado por: \[ \bf\tilde{S}=\bf\tilde{S}_1+\bf\tilde{S}_2=\frac{\bf{S}_1}{n_1}+\frac{\bf{S}_2}{n_2} \] donde \(\bf S_i\) es la matrix de varianzas y covarianzas para \(i=1,\,2\).

Bajo la hipótesis nula se tiene que:

\[T^2\sim\,q\,F_{p,\nu}\] Donde \[q=p+2D-\frac{6D}{p\left(p-1\right)+2}\] \[\nu=\frac{p\left(p+2\right)}{3D}\]

\[ D=\frac{1}{2} \sum^2_{i=1} \frac{1}{n_i-1} \left\lbrace tr\left[\left(\bf I-\left(\bf\tilde{S}^{-1}_1+\bf\tilde{S}^{-1}_2\right)^{-1}\bf\tilde{S}^{-1}_i\right)^2\right]+ \left[tr\left(\bf I-\left(\bf\tilde{S}^{-1}_1+\bf\tilde{S}^{-1}_2\right)^{-1}\bf\tilde{S}^{-1}_i\right)\right]^2 \right\rbrace \]

Generación de las muestras aleatorias

Para generar muestras aleatorias desde una distribución normal bivariada puede ser usada la siguiente función:

gen_data <- function(mu1, sigma1, n1, mu2, sigma2, n2) {
  dt1 <- MASS::mvrnorm(n1, mu=mu1, Sigma=sigma1)
  dt1 <- cbind(dt1, 1)
  dt2 <- MASS::mvrnorm(n2, mu=mu2, Sigma=sigma2) 
  dt2 <- cbind(dt2, 2)
  datos <- as.data.frame(rbind(dt1, dt2))
  colnames(datos) <- c('x', 'y', 'pop')
  datos$pop <- as.factor(datos$pop)
  return(datos)
}

Por ejemplo, para generar dos muestras aleatorias normales bivariadas se utiliza:

n1 <- 45
mu1 <- c(200, 500)
sigma1 <- matrix(c(13825.3, 23823.4,
               23823.4, 73107.4), ncol=2)

n2 <- 55
mu2 <- c(300, 700)
sigma2 <- matrix(c(8632.0, 19616.7,
               19616.7, 55964.5), ncol=2)

datos <- gen_data(mu1, sigma1, n1, mu2, sigma2, n2)

library(ggplot2)
ggplot(data=datos, aes(x=x, y=y, color=pop)) +
  geom_point() +
  annotate("point", x=mu1[1], y=mu1[2], shape=8, size=5) +
  annotate("point", x=mu2[1], y=mu2[2], shape=8, size=5)

Estudio de simulación

Asumiendo diferentes distancias (\(\sqrt{2}\delta\)) entre los vectores de medias de las dos poblaciones se construyen las curvas de potencia para la Prueba de Johansen con diferentes tamaños de muestra (\(n\)). Para generar las muestras normales bivariadas para \(\delta\) y \(n\) específicos se utiliza la función:

auxiliar <- function(delta, n) {
  mu1 <- c(200, 500)
  sigma1 <- matrix(c(13825.3, 23823.4,
                   23823.4, 73107.4), ncol=2)
  mu2 <- mu1 + delta
  sigma2 <- matrix(c(8632.0, 19616.7,
                   19616.7, 55964.5), ncol=2)
  dat <- gen_data(mu1, sigma1, n1=n, mu2, sigma2, n2=n)
  return(dat)
}

El P-valor de la prueba de Johansen para las muestra generadas está dado por:

one_simul <- function(delta, n) {
  datos <- auxiliar(delta, n)
  medias <- lapply(split(x=datos[, 1:2], f=datos$pop, drop=TRUE), colMeans)
  varia <- lapply(split(x=datos[, 1:2], f=datos$pop, drop=TRUE), var)
  N <- lapply(split(x=datos[, 1:2], f=datos$pop, drop=TRUE), nrow)
  res <- two_mean_vector_test(xbar1=medias$`1`, s1=varia$`1`, n1=N$`1`, 
                              xbar2=medias$`2`, s2=varia$`2`, n2=N$`2`,
                              method = 'johansen')
  res$p.value
}

Para estimar la potencia de la prueba para un \(\delta\) y \(n\) específico se utiliza:

mult_simul <- function(x) {
  delta <- x[1]
  n <- x[2]
  nrep <- x[3]
  res <- replicate(n=nrep, expr=one_simul(delta=delta, n=n))
  mean(res < 0.05)
}

Los valores de \(\delta\), \(n\) y el número de repeticiones a considerar son:

delta <- seq(from=0, to=150, by=10)
n <- seq(from=10, to=150, by=10)
nrep <- 5000
casos <- expand.grid(delta=delta, n=n, nrep=nrep)

La tasa de rechazo para los diferentes casos considerados se obtiene a partir de:

tasa <- apply(X=casos, MARGIN=1, FUN=mult_simul)
res <- cbind(casos, tasa)
res$size <- as.factor(res$n)

Un gráfico con los resultados es:

ggplot(data=res, aes(x=delta, y=tasa, group=size)) +
  geom_line(aes(color=size)) + 
  labs(x=expression(delta), y="Potencia")