Hipótesis a estudiar

H01: Los dos perfiles son similares (paralelos)

\[\forall k=1...p-1,\hspace{0.3cm}\mu_{tnt,k+1}-\mu_{tnt,k}=\mu_{c,k+1}-\mu_{c,k}\]

H02: Los dos perfiles están al mismo nivel \[\frac{1}{p}\sum_{j=1}^{p}{\mu_{tnt,j}}=\frac{1}{p}\sum_{j=1}^{p}{\mu_{c,j}}\]

H03: No hay diferencia entre las medias de los tests (no hay efecto principal por test) \[\frac{\mu_{tnt,1}+\mu_{c,1}}{2}=\frac{\mu_{tnt,2}+\mu_{c,2}}{2}=...=\frac{\mu_{tnt,p}+\mu_{c,p}}{2}\]

El gráfico de los perfiles es:

Resumen del Ejercicio P3.1 - 4

No se había rechazado H01, con un p-valor de 0.5673. Entonces se había estudiado la hipótesis H03 que se había rechazado con un p-valor de \(\approx\) 0

Ejercicio P3.2 - 10 Punto a

¿es necesario suponer que las matrices de covarianza en las dos poblaciones sean las mismas o se puede hacer un supuesto más débil?

Para el test de igualdad de medias entre dos poblaciones, el supuesto de igualdad de las matrices de covarianza se puede relajar sí:

  1. si las dos poblaciones son del mismo tamaño:
    • \(n_1=n_2=n\) y \(n\) suficientemente grande, entonces el estadístico \(T_0^2\) es asymptoticamente insensible a \(\Sigma_1 \neq \Sigma_2\).
    • se puede correr un test sobre la media de la variable \(y=x_1-x_2\) que tiene una distribución \(N_p(\mu_1-\mu_2,\Sigma_1 + \Sigma_2)\). Se testea \(H_0 : \mu_y=0\). El estadístico resultante tiene distribución \(T_{p,n-1}^2\) mientras que en la solución previa, el estadístico se basa en \(2(n-1)\) grados de libertad.
  2. si \(n_1\approx n_2\) y son muy grandes, \(\Sigma_1 \neq \Sigma_2\) no impacta el nivel de significación y el poder del test.

Sino hay que usar una solución al problema de Fisher-Behrens, con la de James, Yao o un bootstrap paramétrico.

En nuestro caso tenemos \(n_1=n_2=n=27\). Si consideramos que 27 es suficientemente grande y buscamos aplicar la solución 1a. Primero notamos que en el caso \(n_1=n_2=n\) y \(\Sigma_1=\Sigma_2\), bajo \(H_0\) la distribución del estadístico de Hotelling para dos muestras tiende a una \(\chi^2_p\):

\(T_0^2= \frac{(n-2)p}{n-p-1}F_{p,n-p-1}\)

cuando \(n\longrightarrow+\infty\)

\[\frac{(n-2)p}{n-p-1}\longrightarrow p\]

\[T_0^2= pF_{p,n-p-1 \approx }\chi^2_p\]

Cuando \(n_1=n_2=n\) es suficientemente grande se puede suponer \(S_1=\Sigma_1\) y \(S_2=\Sigma_2\). Entonces bajo, \(H_0\) el estadístico de de Hotelling para dos muestras:

\[T_0^2=\frac{n_1n_2}{n_1+n_2}(\bar{x_1}-\bar{x_2})^T(\frac{(n_1-1)S_1+(n_2-2)S_2}{n_1+n_2-2})^{-1}(\bar{x_1}-\bar{x_2})\]

\[T_0^2=\frac{n_1+n_2-2}{n_1+n_2}(\bar{x_1}-\bar{x_2})^T(\frac{n_1-1}{n_1n_2}S_1+\frac{n_2-1}{n_1n_2}{S_2})^{-1}(\bar{x_1}-\bar{x_2})\]

\[T_0^2=(\bar{x_1}-\bar{x_2})^{T}(\frac{1}{n_2} \Sigma_1+\frac{1}{n_1} \Sigma_2)^{-1}(\bar{x_1}-\bar{x_2})\]

Pero:

\(\bar{x_1}-\bar{x_2} \sim N_p(0,\frac{1}{n_1}\Sigma_1+\frac{1}{n_2}\Sigma_2)\)

Entonces:

\((\bar{x_1}-\bar{x_2})^T(\frac{1}{n_1}\Sigma_1+\frac{1}{n_2}\Sigma_2)^{-1}(\bar{x_1}-\bar{x_2}) \sim \chi^2_p\)

Cuando \(n_1=n_2=n\) suficientemente grande, el estadístico T tiene la distribución asintótica correcta.

Exprese matemáticamente la hipótesis sobre la igualdad de covarianzas que sea necesaria en cada caso

En los tests para H01, H02 y H03, se testea la igualdad de medias de transformaciones de \(x_1\) y \(x_2\).

Para H01:
Definimos \(C_1\):

\[ C_1=\left(\begin{array}{cccccccccc} 1&-1&0&0&0&0&0&0&0&0\\ 0&1&-1&0&0&0&0&0&0&0\\ 0&0&1&-1&0&0&0&0&0&0\\ 0&0&0&1&-1&0&0&0&0&0\\ 0&0&0&0&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&0&0&0\\ 0&0&0&0&0&0&1&-1&0&0\\ 0&0&0&0&0&0&0&1&-1&0\\ 0&0&0&0&0&0&0&0&1&-1 \end{array}\right) \]

\(y_1=C_1x_1, \hspace{0.5cm}y_1\sim N_p(C_1\mu_1,C_1\Sigma_1C_1^T)\)

\(y_2=C_1x_2, \hspace{0.5cm}y_2\sim N_p(C_1\mu_2,C_1\Sigma_2C_1^T)\)

H01 equivale es testear: \(\mu_{y1}=\mu_{y2}\)

La hipótesis sobre covarianzas necesaria es entonces \[\Sigma_{y1}=\Sigma_{y2}\iff C_1\Sigma_1C_1^T=C_1\Sigma_2C_1^T\]

Para H02:
Definimos \(C_2\):

\[ \begin{align} C_2 &= (1,1,1,1,1,1,1,1,1,1) \end{align} \]

\(y_1=C_2x_1, \hspace{0.5cm}y_1\sim N_p(C_2\mu_1,C_2\Sigma_1C_2^T)\)

\(y_2=C_2x_2, \hspace{0.5cm}y_2\sim N_p(C_2\mu_2,C_2\Sigma_2C_2^T)\)

H02 equivale es testear: \(\mu_{y1}=\mu_{y2}\)

La hipótesis sobre covarianzas necesaria es entonces \[\Sigma_{y1}=\Sigma_{y2}\iff C_2\Sigma_1C_2^T=C_2\Sigma_2C_2^T\iff\sum_{i,j=1}^{p}{\sigma_{1,i,j}}=\sum_{i,j=1}^{p}{\sigma_{2,i,j}}\]

Para H03:
Definimos \(C_1\): \[ C_1=\left(\begin{array}{cccccccccc} 1&-1&0&0&0&0&0&0&0&0\\ 0&1&-1&0&0&0&0&0&0&0\\ 0&0&1&-1&0&0&0&0&0&0\\ 0&0&0&1&-1&0&0&0&0&0\\ 0&0&0&0&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&0&0&0\\ 0&0&0&0&0&0&1&-1&0&0\\ 0&0&0&0&0&0&0&1&-1&0\\ 0&0&0&0&0&0&0&0&1&-1 \end{array}\right) \]

\(y_1=C_1x_1, \hspace{0.5cm}y_1\sim N_p(C_1\mu_1,C_1\Sigma_1C_1^T)\)

\(y_2=-C_1x_2, \hspace{0.5cm}y_2\sim N_p(-C_1\mu_2,C_1\Sigma_2C_1^T)\)

H03 equivale es testear: \(\mu_{y1}=\mu_{y2}\)

La hipótesis sobre covarianzas necesaria es entonces \[\Sigma_{y1}=\Sigma_{y2}\iff C_1\Sigma_1C_1^T=C_1\Sigma_2C_1^T\]

Ejercicio P3.2 - 10 Punto b

Para cada una de las situaciones de interés estudie si el supuesto de igualdad de covarianzas es razonable.

Para H01 y H03
Para facilitar el cáclulo en R usamos que como \(n_1=n_2=n\):

\[\gamma^*=\frac{|\frac{Q_1}{n_1}|^{\frac{n_1}{2}}|\frac{Q_2}{n_2}|^{\frac{n_2}{2}}}{|\frac{U}{n_1+n_2}|^{\frac{n_1+n_2}{2}}}=\frac{|\frac{Q_1}{n}|^{\frac{n}{2}}|\frac{Q_2}{n}|^{\frac{n}{2}}}{|\frac{U}{2n}|^{\frac{2n}{2}}}=\left(\frac{|\frac{Q_1}{n}||\frac{Q_2}{n}|}{|\frac{U}{2n}|}\right)^{\frac{n}{2}}\frac{1}{|\frac{U}{2n}|^{\frac{n}{2}}}\]

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/(2*n1)))^(n1/2))/(det(U/(2*n1))^(n1/2)))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H01<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)

Obtenemos un p-valor de 0.2119, entonces no rechazamos al 1%.

Para H02

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rep(1,10)
sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/(2*n1)))^(n1/2))/(det(U/(2*n1))^(n1/2)))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H02<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)

Obtenemos un p-valor de 0.2455, entonces no rechazamos al 1%.

Verificamos el supuesto de normalidad

Corriendo un test de Shapiro multivariado sobre la muestra bruta, obtenemos un p-valor de 0.0305 para la muestra C. Obtenemos un p-valor de 0.1981 para la muestra TNT. Rechazo la normalidad al 5% para la muestra C y rechazo al 20% para la muestra TNT.

En realidad lo que importa es la normalidad de las variables transformadas que son las que testeamos efectivamente. Pasamos entonces a testear normalidad para las variables transformadas.

Para H01 y H03

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
library(mvShapiroTest)
s_res1<-mvShapiro.Test(as.matrix(sample1)) 
s_res2<-mvShapiro.Test(as.matrix(sample2))

Corriendo un test de Shapiro multivariado, obtenemos un p-valor de 0.4624 para la muestra C. Obtenemos un p-valor de 0.9707 para la muestra TNT. En ambos casos no rechazo la normalidad al 20%.

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)

C<-rep(1,10)
sample1<-(t(C%*%t(sample1)))[,1]
sample2<-(t(C%*%t(sample2)))[,1]

#Puedo correr un test univariado porque son vectores
s_res1<-shapiro.test(sample1)
s_res2<-shapiro.test(sample2)

Corriendo un test de Shapiro multivariado, obtenemos un p-valor de 0.3746 para la muestra C. Obtenemos un p-valor de 0.8619 para la muestra TNT. En ambos casos no rechazo la normalidad al 20%.

Era válido aplicar el test de igualdad de matrices de covarianzas para H01, H02 y H03.

TODO LOS QUE ESTÁ ABAJO NO SIRVE

No sirve porque no rechazamos normalidad para las variables que testeamos efectivamente. Lo dejo por si el código les sirve.

Calculamos un estimador robusto de la matriz de covarianza de la muestra C usando covRob de RobStatTM. Realizamos un QQ-plot de las distancias de Mahalanobis robustas obtenidas contra una distribución Chi-cuadrado con 10 grados de libertad.

Quitamos los outliers identificados. Volvemos a graficar QQ-Plot y Boxplot hasta obtener los gráficos siguientes donde quedan un sólo outlier en cada muestra

Corriendo un test de Shapiro multivariado de nuevo, obtenemos un p-valor de 0.5977 para la muestra C. Obtenemos un p-valor de 0.4469 para la muestra TNT. No rechazamos normalidad al 20%

Vuelvo a correr los tests sobre covarianzas sin los outliers.

Para H01 y H03

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample1<-sample1[-c(18,19,26),]
sample1<-sample1[-c(10,11,7,13),]
sample1<-sample1[-c(2),]
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
sample2<-sample2[-c(5,6,15,20,22,23),]
sample2<-sample2[-c(18),]
sample2<-sample2[-c(9),]

C<-rbind(c(1,-1,0,0,0,0,0,0,0,0),
         c(0,1,-1,0,0,0,0,0,0,0),
         c(0,0,1,-1,0,0,0,0,0,0),
         c(0,0,0,1,-1,0,0,0,0,0),
         c(0,0,0,0,1,-1,0,0,0,0),
         c(0,0,0,0,0,1,-1,0,0,0),
         c(0,0,0,0,0,0,1,-1,0,0),
         c(0,0,0,0,0,0,0,1,-1,0),
         c(0,0,0,0,0,0,0,0,1,-1)
)

sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

nm<-min(n1,n2)
n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/n))^(nm/2))*((det(Q_1/n1)^((n1-nm)/2))*(det(Q_2/n2)^((n2-nm)/2))/(det(U/n)^((n-nm)/2))))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)
p_valor_H01<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
p_valor_H01
## [1] 8.930559e-08

Obtenemos un p-valor de approx 0, entonces rechazamos al 1%.

Para H02

remove(list=ls())
sample1<-read.table("ej4pra3-c.txt",row.names = 1)
sample1<-sample1[-c(18,19,26),]
sample1<-sample1[-c(10,11,7,13),]
sample1<-sample1[-c(2),]
sample2<-read.table("ej4pra3-tnt.txt",row.names = 1)
sample2<-sample2[-c(5,6,15,20,22,23),]
sample2<-sample2[-c(18),]
sample2<-sample2[-c(9),]

C<-rep(1,10)
sample1<-t(C%*%t(sample1))
sample2<-t(C%*%t(sample2))
# Test hipotesisis de covarianzas para H01
n1<-nrow(sample1)
n2<-nrow(sample2)
x_barra1<-colMeans(sample1)
x_barra2<-colMeans(sample2)
p<-qr(sample1)$rank

Q_1<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample1)){
  Q_1<-Q_1+(sample1[j,]-x_barra1)%*%t(sample1[j,]-x_barra1)
}


Q_2<-matrix(rep(0,p*p),ncol=p)
for (j  in 1:nrow(sample2)){
  Q_2<-Q_2+(sample2[j,]-x_barra2)%*%t(sample2[j,]-x_barra2)
}


U<-Q_1+Q_2

k<-2

nm<-min(n1,n2)
n<-n1+n2
estadis<-log(((det(Q_1/n1)*det(Q_2/n2)/det(U/n))^(nm/2))*((det(Q_1/n1)^((n1-nm)/2))*(det(Q_2/n2)^((n2-nm)/2))/(det(U/n)^((n-nm)/2))))*(-2)

nu<-k*p*(p+1)/2
nu_1<-p*(p+1)/2

chi_alpha<-qchisq(0.01,nu-nu_1,lower.tail = FALSE)

p_valor_H02<-pchisq(estadis,nu-nu_1,lower.tail = FALSE)
p_valor_H02
## [1] 0.4693055

Obtenemos un p-valor de approx.0.4693, entonces no rechazamos al 1%.