Ejercicio 2:


a) Se realiza el cálculo del estadístico \(T^2\) de Hotelling para testear la igualdad entre los vectores de medias de ambos grupos, basado en \(y = (y_1; y_2)^t\). Se verifica que se obtiene el mismo valor que el del estadístico de Hotelling si se usan las variables originales. (Comparando con el ejercicio 3 de la práctica 2). Se concluye que se rechaza la hipótesis nula. 


library(readxl)
data<-read_xlsx("D:/Maestria en Estadistica 2015/AnalisisMultivariadoI2019/PracticasCorregidas/Datos_Ej2Prac3_2017.xlsx",col_names = FALSE)
data<-data.frame(data)
colnames(data)<-c("Tipo","Antenas","Alas")

x1<-as.numeric(data$Antenas)
x2<-as.numeric(data$Alas)

y1<-x1+x2
y2<-x2

X<-cbind(x1,x2)
Y<-cbind(y1,y2)
Sx<-cov(X)
Sy<-cov(Y)


mu1_y<-apply(Y[1:9,],2,mean)
mu2_y<-apply(Y[10:15,],2,mean)
mu1_x<-apply(X[1:9,],2,mean)
mu2_x<-apply(X[10:15,],2,mean)

n<-dim(Y)[1]
p<-dim(Y)[2]
n1<-sum(data$Tipo=="Af")
n2<-n-n1

S1<-cov(Y[1:9,])
S2<-cov(Y[10:15,])
S<-((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)

T0<-n1*n2/(n1+n2)*t(mu1_y-mu2_y)%*% solve(S)%*%(mu1_y-mu2_y)#55.8807
F0<-T0*(n1+n2-p-1)/(p*(n1+n2-2))#25.79109
falpha = qf(0.05, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#3.885294
pvalor = pf(F0, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#4.519337e-05 (RECHAZO)

falpha*p*(n1+n2-2)/(n1+n2-2-2+1)#8.418137 < T0 (RECHAZO)
## [1] 8.418137
#Lo pruebo con x1 y x2 para ver que da igual que antes con y1 e y2:

S1_x<-cov(X[1:9,])
S2_x<-cov(X[10:15,])
S_x<-((n1-1)*S1_x+(n2-1)*S2_x)/(n1+n2-2)
T0<-n1*n2/(n1+n2)*t(mu1_x-mu2_x)%*% solve(S_x)%*%(mu1_x-mu2_x)#55.8807
F0<-T0*(n1+n2-p-1)/(p*(n1+n2-2))#25.79109
falpha = qf(0.05, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#3.885294
pvalor = pf(F0, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#4.519337e-05 (RECHAZO)

falpha*p*(n1+n2-2)/(n1+n2-2-2+1)#8.418137 < T0 (RECHAZO)
## [1] 8.418137
#pruebo con alpha=0.0025
falpha = qf(0.0025, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#10.28651
pvalor = pf(F0, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#4.519337e-05 (RECHAZO)

falpha*p*(n1+n2-2)/(n1+n2-2-2+1)#22.28743 < T0 (RECHAZO)
## [1] 22.28743


b) En los tests de t univariados de nivel de significación 0,05 realizados sobre cada variable por separado se aceptaría la hipótesis de igualdad de medias. 


#test univariado- Variable longitud de antenas + longitud de alas
S1.y<-sum((Y[1:9,1]-mu1_y[1])^2)/(n1-1)
S2.y<-sum((Y[10:15,1]-mu2_y[1])^2)/(n2-1)
S.1.y<-((n1-1)*S1.y+(n2-1)*S2.y)/(n1+n2-2)
t_obs<-abs((mu1_y[1]-mu2_y[1])/(sqrt(S.1.y)*sqrt(1/n1+1/n2)))#0.6609713 
t_obs
##        y1 
## 0.6609713
t<-qt(0.025,n1+n2-2, ncp=0, lower.tail = FALSE, log.p = FALSE)#2.160369
t
## [1] 2.160369
#como t_obs<t NO Rechazo!
pvalor = pt(t_obs, n1+n2-2, ncp=0, lower.tail = FALSE, log.p = FALSE)
pvalor#0.2600865 (NO Rechazo!)
##        y1 
## 0.2600865
#test univariado- Variable longitud de alas
S1.y<-sum((Y[1:9,2]-mu1_y[2])^2)/(n1-1)
S2.y<-sum((Y[10:15,2]-mu2_y[2])^2)/(n2-1)
S.2.y<-((n1-1)*S1.y+(n2-1)*S2.y)/(n1+n2-2)
t_obs<-abs((mu1_y[2]-mu2_y[2])/(sqrt(S.2.y)*sqrt(1/n1+1/n2)))#2.004721 
t_obs
##       y2 
## 2.004721
t<-qt(0.025,n1+n2-2, ncp=0, lower.tail = FALSE, log.p = FALSE)#2.160369
t
## [1] 2.160369
#como t_obs<t NO Rechazo!
pvalor = pt(t_obs, n1+n2-2, ncp=0, lower.tail = FALSE, log.p = FALSE)
pvalor#0.03313874>0.025 (NO Rechazo!)
##         y2 
## 0.03313874


c) Scatterplot de \(y_1\) versus \(y_2\). 


par(pty="s")
plot(y1,y2,pch=c(rep(16,9),rep(15,6)),col=c(rep(1,9),rep(2,6)),xlim=c(2.7,3.8),ylim=c(1.55,2.2),asp=1)
points(mu1_y[1],mu1_y[2],col=1,pch="*")
points(mu2_y[1],mu2_y[2],col=2,pch="*")

ProyPsobreL<-function(P,v,Q){
  R<-as.numeric((P-Q)%*%v/sum(v^2))*v+Q
  return(R)
}

a_hat<-solve(S)%*%(mu1_y-mu2_y)

Proy<-matrix(NA,nrow=n,ncol=p)
for (k in 1:n){
  Proy[k,]<-ProyPsobreL(Y[k,],a_hat,c(3,1.5))
}

r1<-c(-100:100)*a_hat[1]+3
r2<-c(-100:100)*a_hat[2]+1.5
lines(r1,r2,col=3)

points(Proy[1:9,],col=4,pch=16)
points(Proy[10:15,],col="hotpink",pch=15)


d) Se observa que la elipse de confianza de nivel 98% para el vector de diferencia de medias no cubre al vector 0. Sin embargo, el rectángulo correspondiente a los intervalos de confianza univariados de nivel 99% (\((0.99)^2\cong~0.98\)) para las diferencias entre las medias de \(y_1\) e \(y_2\) sí lo contienen. 


A<-solve(S)

eig<-eigen(A)
autovalores<-eig$values
U<-eig$vectors
N<-n1+n2-2
falpha0.02 = qf(0.02, p, n1+n2-p-1, ncp=0, lower.tail = FALSE, log.p = FALSE)#5.516299
c_cuad<-(n1+n2)/(n1*n2)*p*N/(N-p+1)*falpha0.02
c<-sqrt(c_cuad)
y.1<-vector()
y.2<-vector()
tita<-seq(0,2*pi,2*pi/100)

for (i in 1:length(tita)){
  y1<-(c/sqrt(autovalores[1]))*cos(tita[i])
  y2<-(c/sqrt(autovalores[2]))*sin(tita[i])
  V<-U%*%c(y1,y2)
  y.1[i]<-V[1]
  y.2[i]<-V[2]
}
b<-mu1_y-mu2_y
y.1<-y.1+b[1]
y.2<-y.2+b[2]
par(pty="s")
plot(y.1,y.2,type="l",lwd=2,col="orange",xlab="y1",ylab="y2")
points(0,0,pch=15)


t_0.005<-qt(0.005,n1+n2-2, ncp=0, lower.tail = FALSE, log.p = FALSE)#3.012276
IC1_a<-(mu1_y[1]-mu2_y[1])-(sqrt(S.1.y)*sqrt(1/n1+1/n2))*t_0.005#-0.2292513
IC1_b<-(mu1_y[1]-mu2_y[1])+(sqrt(S.1.y)*sqrt(1/n1+1/n2))*t_0.005#0.3581401


IC2_a<-(mu1_y[2]-mu2_y[2])-(sqrt(S.2.y)*sqrt(1/n1+1/n2))*t_0.005#-0.3058722
IC2_b<-(mu1_y[2]-mu2_y[2])+(sqrt(S.2.y)*sqrt(1/n1+1/n2))*t_0.005#0.06142779

segments(IC1_a,IC2_a,IC1_b,IC2_a,lty=2,col=4,lwd=2)
segments(IC1_a,IC2_b,IC1_b,IC2_b,lty=2,col=4,lwd=2)
segments(IC1_a,IC2_a,IC1_a,IC2_b,lty=2,col=6,lwd=2)
segments(IC1_b,IC2_a,IC1_b,IC2_b,lty=2,col=6,lwd=2)
abline(h=0,lty=2)
abline(v=0,lty=2)