#1.- INTRODUCCIÓN

A continuación se presenta una serie de ejercicios propuesto que permitiran interiorizar el análisis de datos multivariados utilizandos las técnicas ANOVA y MANOVA.

library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
library(mvShapiroTest)
library(biotools)
## Warning: package 'biotools' was built under R version 4.2.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.2.3
## ---
## biotools version 4.2
library(outliers)
library(ICSNP)
## Warning: package 'ICSNP' was built under R version 4.2.3
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 4.2.3
## Loading required package: ICS
## Warning: package 'ICS' was built under R version 4.2.3
library(mvtnorm)
library(ICS)
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

#2.- DATA SET

Los ejercicios mostrados son desarrollados a partir de los siguientes datos correspondientes a una toma realizada con un espectrorradiómetro en los anchos de bandas correspondientes a 560 nm y 720 nm de las especies SS, JL y LP.

library(readxl)
Datos_ANOVA_MANOVA <- read_excel("C:/Users/Usuario/Desktop/Luis/Data Science/REPASO/R/Clase 14/Datos_ANOVA_MANOVA.xlsx")
View(Datos_ANOVA_MANOVA)
Reflectancia<-Datos_ANOVA_MANOVA
head(Reflectancia)
## # A tibble: 6 × 3
##   D_560nm D_720nm Species
##     <dbl>   <dbl> <chr>  
## 1    9.33    19.1 SS     
## 2    8.74    19.6 SS     
## 3    9.31    19.2 SS     
## 4    8.27    16.4 SS     
## 5   10.2     25   SS     
## 6   10.1     25.3 SS

#3.- MANOVA

ggplot(data=Reflectancia, aes(Species,D_560nm))+
  geom_point()+ggtitle("Comportamiento Firma 720 nm") + xlab("Especies")+ylab("Reflectancia")

ggplot(data=Reflectancia, aes(D_560nm, D_720nm, color=Species))+
  geom_point()+ggtitle("Comportamiento Firma 560 nm va 720 nm") 

MN = manova(cbind(Reflectancia$D_560nm, Reflectancia$D_720nm)~Reflectancia$Species)
MN
## Call:
##    manova(cbind(Reflectancia$D_560nm, Reflectancia$D_720nm) ~ Reflectancia$Species)
## 
## Terms:
##                 Reflectancia$Species Residuals
## resp 1                       965.181  2147.714
## resp 2                      2026.856  7536.997
## Deg. of Freedom                    2        33
## 
## Residual standard errors: 8.067357 15.11271
## Estimated effects may be unbalanced
summary(MN, test="Wilks")
##                      Df   Wilks approx F num Df den Df Pr(>F)  
## Reflectancia$Species  2 0.67704   3.4452      4     64  0.013 *
## Residuals            33                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obsérvese que se presentan diferencias significativas entre las dos especies.

#4.- ANOVAS

ANOVA1 <- aov(Reflectancia$D_560nm~Reflectancia$Species)
summary(ANOVA1)
##                      Df Sum Sq Mean Sq F value  Pr(>F)   
## Reflectancia$Species  2  965.2   482.6   7.415 0.00219 **
## Residuals            33 2147.7    65.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ANOVA1)

El estadístico F asciende a 7.41 con un p valor de 0.002 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Hay diferencias significtivas por especies.

res1 <- ANOVA1$residuals
shapiro.test(res1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.88653, p-value = 0.001491

No se acepta la hipótesis de normalidad.

ANOVA2 <- aov(Reflectancia$D_720nm~Reflectancia$Species)
summary(ANOVA2)
##                      Df Sum Sq Mean Sq F value Pr(>F)  
## Reflectancia$Species  2   2027  1013.4   4.437 0.0196 *
## Residuals            33   7537   228.4                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(ANOVA2)

El estadístico F asciende a 4.43 con un p valor de 0.019 < 0.05. Tenemos evidencia empírica suficiente para rechazar H0 al 5%. Hay diferencias significtivas por especies.

res2<-ANOVA2$residuals
shapiro.test(res2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.81542, p-value = 3.347e-05

No se acepta la hipótesis de normalidad.

#5.- Correlación

cor(Reflectancia$D_560nm, Reflectancia$D_720nm)
## [1] 0.8130816

La correlación es fuerte y positiva.

cor.test(Reflectancia$D_560nm, Reflectancia$D_720nm)
## 
##  Pearson's product-moment correlation
## 
## data:  Reflectancia$D_560nm and Reflectancia$D_720nm
## t = 8.144, df = 34, p-value = 1.691e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6611588 0.9009499
## sample estimates:
##       cor 
## 0.8130816

El p valor es 0.000 < 0.05. Por tanto, la correlación es estadísticamente significativa.

#6.- Normalidad

shapiro.test(Reflectancia$D_560nm)
## 
##  Shapiro-Wilk normality test
## 
## data:  Reflectancia$D_560nm
## W = 0.64673, p-value = 4.538e-08

No se acepta la normalidad.

shapiro.test(Reflectancia$D_720nm)
## 
##  Shapiro-Wilk normality test
## 
## data:  Reflectancia$D_720nm
## W = 0.84251, p-value = 0.0001286

No se acepta normalidad.

library(mvShapiroTest)
mvShapiro.Test(as.matrix(Reflectancia[,1:2]))
## 
##  Generalized Shapiro-Wilk test for Multivariate Normality by
##  Villasenor-Alva and Gonzalez-Estrada
## 
## data:  as.matrix(Reflectancia[, 1:2])
## MVW = 0.81281, p-value = 3.948e-08

No se acepta la normalidad multivariada.

#7Homogeneidad de varianzas o esfericidad

bartlett.test(Reflectancia$D_560nm~Reflectancia$Species)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Reflectancia$D_560nm by Reflectancia$Species
## Bartlett's K-squared = 32.714, df = 2, p-value = 7.876e-08

Hay heterocedasticidad.

bartlett.test(Reflectancia$D_720nm~Reflectancia$Species)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Reflectancia$D_720nm by Reflectancia$Species
## Bartlett's K-squared = 1.7462, df = 2, p-value = 0.4177

Tenemos homocedasticidad, dado que p valor es 0.4177 > 0.05.

library(biotools)
boxM(as.matrix(Reflectancia[,-3]), as.matrix(Reflectancia$Species))
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  as.matrix(Reflectancia[, -3])
## Chi-Sq (approx.) = 41.969, df = 6, p-value = 1.865e-07

Las matrices de varianzas y covarianzas no son similares.

#8.- Atípicos

ggplot(data=Reflectancia, aes(Species,D_560nm))+ geom_boxplot()+ggtitle("Comportamiento Firma 560 nm") + xlab("Especies")+ylab("Reflectancia")

ggplot(data=Reflectancia, aes(Species,D_720nm))+ geom_boxplot()+ggtitle("Comportamiento Firma 720 nm") + xlab("Especies")+ylab("Reflectancia")

dixon.test(sample(Reflectancia$D_560nm, size = 30))
## 
##  Dixon test for outliers
## 
## data:  sample(Reflectancia$D_560nm, size = 30)
## Q = 0.20418, p-value = 0.8265
## alternative hypothesis: highest value 44.74 is an outlier

El p valor es mayor que 0.05. Los atípicos no son influyentes.

dixon.test(sample(Reflectancia$D_720nm, size = 30))
## 
##  Dixon test for outliers
## 
## data:  sample(Reflectancia$D_720nm, size = 30)
## Q = 0.56582, p-value < 2.2e-16
## alternative hypothesis: highest value 78.57 is an outlier

El p valor es mayor que 0.05. Los atípicos no son influyentes.

vec.medias=apply(Reflectancia[,1:2],2,mean);vec.medias
##  D_560nm  D_720nm 
## 14.30139 35.78806
T2=c()
for(j in 1:dim(Reflectancia)[1]){
  T2[j]=c((t(t(Reflectancia[j,1:2])-(vec.medias)))%*%solve(var(Reflectancia[,1:2]))%*%as.matrix(t(Reflectancia[j,1:2])-(vec.medias)))
}
T2
##  [1]  1.26541020  1.09383567  1.24128875  1.67375718  0.45417223  0.41656430
##  [7]  0.27561128  0.34847917  0.04318851  0.07844282  0.10339451 12.43210949
## [13]  0.05790637  0.13704612  0.07555811  0.05807567  0.09033585  0.25119830
## [19]  0.47601009  0.19901130  7.16415107 10.42111846  5.77865326 11.83215879
## [25]  0.57530618  0.82479717  0.68473177  0.71821599  0.39872455  0.69415980
## [31]  0.40375328  0.53175909  0.88021407  0.12249645  0.14489054  8.05347359
LS=qchisq(0.95,df=dim(Reflectancia)-1)
colores=ifelse(T2>LS,"darkred","darkgreen")
plot(T2,col=colores,pch=19,cex=0.85,xlab="Observación")
grid(20,20,col="lightblue")
abline(h=LS)
etiquetas=which(T2>LS)
text(etiquetas,c(LS)+0.4,"outlier")

plot(Reflectancia$D_560nm,Reflectancia$D_720nm,pch=19,col=colores)
grid(20,20,col="lightblue")

#9.- Comparaciones múltiples

# ESPECIE SS - JL

data.a <- filter(Reflectancia, Species != "LP")
PruebaH21 <- HotellingsT2(as.matrix(data.a[,1:2])~as.matrix(data.a[,3]), mu = c(0,0)) ##
PruebaH21
## 
##  Hotelling's two sample T2-test
## 
## data:  as.matrix(data.a[, 1:2]) by as.matrix(data.a[, 3])
## T.2 = 3.2938, df1 = 2, df2 = 21, p-value = 0.05699
## alternative hypothesis: true location difference is not equal to c(0,0)

Se acepta la hipótesis de igualdad de medias entre SS y JL.

# ESPECIE SS - LP

data.b <- filter(Reflectancia, Species != "JL")
PruebaH22 <- HotellingsT2(as.matrix(data.b[,1:2])~as.matrix(data.b[,3]), mu = c(0,0)) ##
PruebaH22
## 
##  Hotelling's two sample T2-test
## 
## data:  as.matrix(data.b[, 1:2]) by as.matrix(data.b[, 3])
## T.2 = 2.4676, df1 = 2, df2 = 21, p-value = 0.109
## alternative hypothesis: true location difference is not equal to c(0,0)

Se acepta la hipótesis de igualdad de medias entre SS y LP.

# ESPECIE JL - LP

data.c <- filter(Reflectancia, Species != "SS")
PruebaH23 <- HotellingsT2(as.matrix(data.c[,1:2])~as.matrix(data.c[,3]), mu = c(0,0)) ##
PruebaH23
## 
##  Hotelling's two sample T2-test
## 
## data:  as.matrix(data.c[, 1:2]) by as.matrix(data.c[, 3])
## T.2 = 4.2569, df1 = 2, df2 = 21, p-value = 0.02806
## alternative hypothesis: true location difference is not equal to c(0,0)

Dado que el p valor es 0.02 < 0.05, en este caso las medias no son iguales entre JL y LP.

# DATOS POR CADA ESPECIE

data.SS <- filter(Reflectancia, Species == "SS")
data.JL <- filter(Reflectancia, Species == "JL")
data.LP <- filter(Reflectancia, Species == "LP")
vec.mediasSS=apply(data.SS[,1:2],2,mean);vec.mediasSS
##  D_560nm  D_720nm 
## 11.54333 30.15167
vec.mediasJL=apply(data.JL[,1:2],2,mean);vec.mediasJL
## D_560nm D_720nm 
## 21.5550 46.3925
vec.mediasLP=apply(data.LP[,1:2],2,mean);vec.mediasLP
##   D_560nm   D_720nm 
##  9.805833 30.820000
varSS=var(data.SS[,1:2]); varSS
##           D_560nm   D_720nm
## D_560nm  9.444588  28.47334
## D_720nm 28.473339 203.87040
varJL=var(data.JL[,1:2]); varJL
##          D_560nm  D_720nm
## D_560nm 178.3878 221.8819
## D_720nm 221.8819 332.1606
varLP=var(data.LP[,1:2]); varLP
##           D_560nm   D_720nm
## D_560nm  7.414354  27.71932
## D_720nm 27.719318 149.15058
# PARÁMETROS DEL INTERVALO

vec.mediasTotal=apply(Reflectancia[,1:2],2,mean);vec.mediasTotal
##  D_560nm  D_720nm 
## 14.30139 35.78806
t1=vec.mediasSS-vec.mediasTotal
t2=vec.mediasJL-vec.mediasTotal
t3=vec.mediasLP-vec.mediasTotal
W=(dim(data.SS)[1]-1)*varSS+(dim(data.JL)[1]-1)*varJL+(dim(data.LP)[1]-1)*varLP; W
##          D_560nm  D_720nm
## D_560nm 2147.714 3058.820
## D_720nm 3058.820 7536.997
pt=qt((0.05/(2*3*(3-1))),df=(dim(Reflectancia)[1]-3)); pt
## [1] -2.806761
t.SS_JL.560=t1[1]-t2[1]
p.SS_JL.560=sqrt(((1/dim(data.SS)[1])+(1/dim(data.JL)[1]))*(W[1,1]/(dim(Reflectancia)[1]-3)))
(limsup.SS_JL.560=t.SS_JL.560+p.SS_JL.560)
##   D_560nm 
## -6.718182
(liminf.SS_JL.560=t.SS_JL.560-p.SS_JL.560)
##   D_560nm 
## -13.30515

IC al 95% está entre -6.71 y -13.30 (-,-). Esto quiere decir que la media SS < media de JL en la variable 560nm.

# INTERVALO ENTRE SS Y JL A UNA LONGITUD DE ONDA DE 720 nm

t.SS_JL.720=t1[2]-t2[2]
p.SS_JL.720=sqrt(((1/dim(data.SS)[1])+(1/dim(data.JL)[1]))*(W[2,2]/(dim(Reflectancia)[1]-3)))
(limsup.SS_JL.720=t.SS_JL.720+p.SS_JL.720)
##  D_720nm 
## -10.0711
(liminf.SS_JL.720=t.SS_JL.720-p.SS_JL.720)
##   D_720nm 
## -22.41057

Los límites son negativos.

IC al 95% entre -10.07 y -22.41 (-,-). Esto quiere decir que la media SS < media JL en la variable 720nm.

# INTERVALO ENTRE SS Y LP A UNA LONGITUD DE ONDA DE 560 nm

t.SS_LP.560=t1[1]-t3[1]
p.SS_LP.560=sqrt(((1/dim(data.SS)[1])+(1/dim(data.LP)[1]))*(W[1,1]/(dim(Reflectancia)[1]-3)))
(limsup.SS_LP.560=t.SS_LP.560+p.SS_LP.560)
##  D_560nm 
## 5.030985
(liminf.SS_LP.560=t.SS_LP.560-p.SS_LP.560)
##   D_560nm 
## -1.555985

El intervalo de confianza es (-,+); por lo tanto, las medias son similares entre SS y LP para la variable 560 nm.

# INTERVALO ENTRE SS Y LP A UNA LONGITUD DE ONDA DE 720 nm

t.SS_LP.720=t1[2]-t3[2]
p.SS_LP.720=sqrt(((1/dim(data.SS)[1])+(1/dim(data.LP)[1]))*(W[2,2]/(dim(Reflectancia)[1]-3)))
(limsup.SS_LP.720=t.SS_LP.720+p.SS_LP.720)
##  D_720nm 
## 5.501403
(liminf.SS_LP.720=t.SS_LP.720-p.SS_LP.720)
##   D_720nm 
## -6.838069

El intervalo de confianza es (-,+); por lo tanto, las medias son similares entre SS y LP para la variable 720 nm.

# INTERVALO ENTRE JL Y LP A UNA LONGITUD DE ONDA DE 560 nm

t.JL_LP.560=t2[1]-t3[1]
p.JL_LP.560=sqrt(((1/dim(data.JL)[1])+(1/dim(data.LP)[1]))*(W[1,1]/(dim(Reflectancia)[1]-3)))
(limsup.JL_LP.560=t.JL_LP.560+p.JL_LP.560)
##  D_560nm 
## 15.04265
(liminf.JL_LP.560=t.JL_LP.560-p.JL_LP.560)
##  D_560nm 
## 8.455682

El IC (+,+), indicando que la media JL > la media de LP en la variable 560 nm.

# INTERVALO ENTRE JL Y LP A UNA LONGITUD DE ONDA DE 720 nm

t.JL_LP.720=t2[2]-t3[2]
p.JL_LP.720=sqrt(((1/dim(data.JL)[1])+(1/dim(data.LP)[1]))*(W[2,2]/(dim(Reflectancia)[1]-3)))
(limsup.JL_LP.720=t.JL_LP.720+p.JL_LP.720)
##  D_720nm 
## 21.74224
(liminf.JL_LP.720=t.JL_LP.720-p.JL_LP.720)
##  D_720nm 
## 9.402764

El IC (+,+), indicando que la media JL > la media de LP en la variable 720 nm.

#10.- CONCLUSIÓN

A partir de los resultados se puede indicar que para las especies SS y JL, JL y LP son diferentes para las dos longitudes de onda mientras que esta característica es contraría para las especies SS y LP, es decir que estos dos sensores no poseen diferencias estadísticas en las dos longitudes de onda.

#*11.- INFORME DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 90 × 3
##    package    loadedversion source        
##    <chr>      <chr>         <chr>         
##  1 biotools   4.2           CRAN (R 4.2.3)
##  2 boot       1.3-28        CRAN (R 4.2.2)
##  3 bslib      0.4.2         CRAN (R 4.2.2)
##  4 cachem     1.0.6         CRAN (R 4.2.2)
##  5 callr      3.7.3         CRAN (R 4.2.3)
##  6 cellranger 1.1.0         CRAN (R 4.2.3)
##  7 cli        3.6.0         CRAN (R 4.2.2)
##  8 colorspace 2.1-0         CRAN (R 4.2.3)
##  9 crayon     1.5.2         CRAN (R 4.2.3)
## 10 DBI        1.1.3         CRAN (R 4.2.3)
## # ℹ 80 more rows