Introducción

Para implementar la prueba de NVM 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 NVM

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 Nel and Van der Merwe (1986). 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\,\frac{\nu p}{\nu-p+1}\,F_{p,\nu-p+1}\] Donde:

\[\nu=\frac{tr\left[\left(\bf\tilde S\right)^2\right]+\left[tr\left(\bf\tilde S\right)\right]^2}{\frac{1}{n_1-1}\left\lbrace tr\left[ \left(\tilde S_1\right)^2\right]+tr\left[ \left(\tilde S_1\right)\right]^2 \right\rbrace+\frac{1}{n_2-1}\left\lbrace tr\left[ \left(\tilde S_2\right)^2\right]+tr\left[ \left(\tilde S_2\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 NVM 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 NVM para las muestras 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 = 'nvm')
  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")