Ejercicio con datos Medidas del caparazón (en milímetros) de las tortugas pintadas
Yeimy Maryury Montilla Montilla; Carlos Eduardo Caceres Gonzalez
02/05/2021
library(readxl)
tortugas <- read_excel("F:/UNAL_MAESTRIA_UNAL/METODOS_MULTIVARIADOS/Trabajo2/datos.xlsx",
sheet = "Carapace_Measurements_ord")
data=data.frame(tortugas)Estadistico T^2
data_Female=subset(data,data$Genero!="Female")
data_Male=subset(data,data$Genero!="Male")
par(mfrow=c(1,2))
boxplot(data_Female[,1:3],ylab='Medida (mm)',main='Hembras', ylim = c(0, 180))
boxplot(data_Male[,1:3],ylab='Medida (mm)',main='Machos', ylim = c(0, 180))
A partir del grafico se puede sugerir que las medidas del caparazón de las tortugas son diferentes entre machos y hembras.
X1_F=tortugas[tortugas$Genero=="Female",1];X1_F=unlist(X1_F)
X1_M=tortugas[tortugas$Genero=="Male",1];X1_M=unlist(X1_M)
X2_F=tortugas[tortugas$Genero=="Female",2];X2_F=unlist(X2_F)
X2_M=tortugas[tortugas$Genero=="Male",2];X2_M=unlist(X2_M)
X3_F=tortugas[tortugas$Genero=="Female",3];X3_F=unlist(X3_F)
X3_M=tortugas[tortugas$Genero=="Male",3];X3_M=unlist(X3_M)dX1=X1_F-X1_M
dX2=X2_F-X2_M
dX3=X3_F-X3_MPrueba \(T^2\) de Hotelling para muestras pareadas.
\[H_0: \begin{bmatrix} \mu_{X1_F} \\\mu_{X2_F} \\ \mu_{X3_F} \end{bmatrix} = \begin{bmatrix} \mu_{X1_M} \\\mu_{X2_M} \\ \mu_{X3_M} \end{bmatrix}\] Diferencias entre variables
df1=data.frame(dX1,dX2,dX3)
df2=t(df1)
colnames(df2)=paste0("obs",1:24)
df2## obs1 obs2 obs3 obs4 obs5 obs6 obs7 obs8 obs9 obs10 obs11 obs12 obs13 obs14
## dX1 5 9 7 4 7 20 19 27 26 21 21 22 22 21
## dX2 7 6 6 2 3 11 12 16 20 13 12 16 8 9
## dX3 1 3 7 3 6 13 7 12 13 11 8 9 8 10
## obs15 obs16 obs17 obs18 obs19 obs20 obs21 obs22 obs23 obs24
## dX1 24 28 29 33 34 30 31 31 31 42
## dX2 14 15 18 14 20 24 19 23 29 26
## dX3 12 16 15 12 21 15 17 18 15 20Vector de media de diferencias
d=rowMeans(df2)
d=as.matrix(d);d## [,1]
## dX1 22.66667
## dX2 14.29167
## dX3 11.33333\[\bar d = \begin{bmatrix} \bar d_{x1} \\ \bar d_{x2} \\ \bar d_{x3} \end{bmatrix} =\begin{bmatrix} 22.67 \\ 14.29 \\ 11.33 \end{bmatrix}\] Matriz de varianzas y covarianzas de las diferencias medias
s_d=cov(df1);s_d #Matriz de varianzas y covarianzas de diff medias## dX1 dX2 dX3
## dX1 101.71014 63.49275 49.07246
## dX2 63.49275 51.78080 32.07246
## dX3 49.07246 32.07246 28.49275\[S_d=\begin{bmatrix} 101.71 & 63.49 & 49.07 \\ 63.49 &51.78 &32.07 \\49.07 & 32.07 &28.49 \end{bmatrix}\]
\(T^2\) Calculado
T2=24*t(d)%*%solve(s_d)%*%d;T2 #calculado## [,1]
## [1,] 122.0209\(T^2\) Tabulado
p=3
n=24
ftabla=qf(p = 0.05,df1 = 3,df2 = 21,lower.tail = F)
TT2=((p*(n-1))/(n-p))*ftabla;TT2 #tabulado## [1] 10.09525Se rechaza la hipotesis nula si \(T^2\) calculado es menor al \(T^2\) tabulado
ifelse(T2>TT2,"Rechazo Ho","No rechazo Ho")## [,1]
## [1,] "Rechazo Ho"De acuerdo con la prueba \(T^2\) y su comparación con \(T^2\) tabulado se rechaza la hipótesis nula. Por tanto, las medidas del caparazón de las tortugas son diferentes según el género.
¿Qué variable influye más en el rechazo de la hipótesis nula?
ICS<-function(d,T2,vari){
li=NULL
ls=NULL
for( i in 1:length(d)){
li[i]=d[i]-sqrt(T2)*sqrt(vari[i,i]/24)
ls[i]=d[i]+sqrt(T2)*sqrt(vari[i,i]/24)
}
return(list(li=li,ls=ls))
}
ICS(d,TT2,s_d)## $li
## [1] 16.125804 9.624673 7.871389
##
## $ls
## [1] 29.20753 18.95866 14.79528\[IC_{x1}=[16.13,29.20] \\IC_{x2}=[9.62,18.96]\\IC_{x3}=[7.87,14.79]\] La variable que más influye en el rechazo de la hipótesis es la variable X1. Por tanto, es la medida donde se presentan mayores diferencias entre machos y hembras.
Manova
Manova = manova(cbind(data$x1,data$x2,data$x3) ~ data$Genero)
summary(Manova, test = 'Wilks')## Df Wilks approx F num Df den Df Pr(>F)
## data$Genero 1 0.38857 23.078 3 44 3.967e-09 ***
## Residuals 46
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De acuerdo al resultado del manova las medidas del caparazón de las tortugas pintadas son diferentes para machos y hembras.
res = Manova$residuals; res## [,1] [,2] [,3]
## 1 -38.04166667 -21.5833333 -14.0416667
## 2 -33.04166667 -18.5833333 -14.0416667
## 3 -33.04166667 -16.5833333 -10.0416667
## 4 -31.04166667 -16.5833333 -10.0416667
## 5 -27.04166667 -14.5833333 -8.0416667
## 6 -13.04166667 -10.5833333 -2.0416667
## 7 -13.04166667 -7.5833333 -6.0416667
## 8 -3.04166667 -3.5833333 -1.0416667
## 9 -3.04166667 -0.5833333 -1.0416667
## 10 -3.04166667 -0.5833333 -1.0416667
## 11 -2.04166667 -2.5833333 -4.0416667
## 12 -0.04166667 -0.5833333 -3.0416667
## 13 1.95833333 -4.5833333 -1.0416667
## 14 1.95833333 -3.5833333 -1.0416667
## 15 4.95833333 2.4166667 0.9583333
## 16 10.95833333 5.4166667 4.9583333
## 17 12.95833333 4.4166667 2.9583333
## 18 16.95833333 4.4166667 3.9583333
## 19 18.95833333 12.4166667 10.9583333
## 20 18.95833333 14.4166667 7.9583333
## 21 21.95833333 12.4166667 9.9583333
## 22 22.95833333 15.4166667 10.9583333
## 23 25.95833333 21.4166667 8.9583333
## 24 40.95833333 29.4166667 14.9583333
## 25 -20.37500000 -14.2916667 -3.7083333
## 26 -19.37500000 -10.2916667 -5.7083333
## 27 -17.37500000 -8.2916667 -5.7083333
## 28 -12.37500000 -4.2916667 -1.7083333
## 29 -11.37500000 -3.2916667 -2.7083333
## 30 -10.37500000 -7.2916667 -3.7083333
## 31 -9.37500000 -5.2916667 -1.7083333
## 32 -7.37500000 -5.2916667 -1.7083333
## 33 -6.37500000 -6.2916667 -2.7083333
## 34 -1.37500000 0.7083333 -0.7083333
## 35 -0.37500000 -0.2916667 -0.7083333
## 36 0.62500000 -2.2916667 -0.7083333
## 37 2.62500000 1.7083333 2.2916667
## 38 3.62500000 1.7083333 0.2916667
## 39 3.62500000 2.7083333 0.2916667
## 40 5.62500000 4.7083333 0.2916667
## 41 6.62500000 0.7083333 -0.7083333
## 42 6.62500000 4.7083333 3.2916667
## 43 7.62500000 6.7083333 1.2916667
## 44 11.62500000 4.7083333 4.2916667
## 45 13.62500000 7.7083333 4.2916667
## 46 14.62500000 6.7083333 4.2916667
## 47 17.62500000 6.7083333 5.2916667
## 48 21.62500000 17.7083333 6.2916667royston::royston.test(res)
## $H
## [1] 0.5829341
##
## $p.value
## [1] 0.5588467
##
## $distribution
## [1] "Data analyzed have a normal distribution."La prueba de normalidad de los residuales indica que siguen una distribución normal.