##Preparación

#Agrego librerias base que se van a necesitar y voy sumando a medida que surjan nuevas

library(readxl)#lectura
library(ggplot2) #
library(dplyr)#manejo de datos
library(tidyverse)#manejo de datos
library(DT)
library(kableExtra)# hacer tablas bonitas
library(lubridate)#manejo de fechas
library(GGally)
library(moments)
library(reshape2)
library(car)
library(plotly)
library(fmsb)
library(geomtextpath)
library(factoextra)
library(flextable)
library(ggrepel)
library(FactoMineR)
library(ggcorrplot)
library(corrplot)

Ejercicio 1

Sea la matriz de varianzas y covarianzas poblacionales:

\[ \begin{equation} \begin{pmatrix} 3 & 1 & 1\\ 1 & 3 & 1\\ 1 & 1 & 5 \end{pmatrix} \end{equation} \]

Correspondiente al vector aleatorio \(X = (X_1,X_2,X_3)\) de media 0.

a) Hallar los autovalores y autovectores de la matriz de varianzas y covarianzas.

mat1 <- matrix(c(3, 1, 1,
             1, 3, 1, 
             1, 1, 5), nrow = 3, byrow = T)

eigen1 <- eigen(mat1)


eigen1[1] 
## $values
## [1] 6 3 2
eigen1[2]$vectors
##            [,1]       [,2]          [,3]
## [1,] -0.4082483 -0.5773503  7.071068e-01
## [2,] -0.4082483 -0.5773503 -7.071068e-01
## [3,] -0.8164966  0.5773503 -1.110223e-16

b) Escribir la expresión de las componentes principales Y = (Y1, Y2, Y3)′ e indique que proporción de la variabilidad explica cada una de ellas.

Si mal no entendí esto es la convinación lineal de las variables con cada uno de los autovectores como coeficientes.

eigen1[2]$vectors[1]*mat1[1,] + eigen1[2]$vectors[2]*mat1[2,] + eigen1[2]$vectors[3]*mat1[3,]
## [1] -2.449490 -2.449490 -4.898979

Validemos usando la funcion prcomp de R:

pca1 <- prcomp(mat1, center=TRUE,scale=TRUE)

pca1
## Standard deviations (1, .., p=3):
## [1] 1.224745e+00 1.224745e+00 4.532467e-17
## 
## Rotation (n x k) = (3 x 3):
##             PC1           PC2        PC3
## [1,] -0.8164966 -2.136625e-17 -0.5773503
## [2,]  0.4082483  7.071068e-01 -0.5773503
## [3,]  0.4082483 -7.071068e-01 -0.5773503

Lo que observamos es que el componente no parece ser el resultado de la combinación lineal. Pero en realidad es LD con el resultado de la función. Y vemos que equivale a \(\lambda*PC_1\)| y esto es correcto ya que \(Av=\lambda v\) que se corresponde a la definición de autovector y autovalor.

NO entiendo porqué tengo que usar scale=TRUE, yo entendía que esto era si usaba la correlacion y no la varianza y covarianza. No entiendo porque dan diferentes y que significa. En particular porque coincide ccuando pongo TRUE, entendía que debía ser al revez

En cuanto a los autovalores deberían darme el porcentaje que explica de los valores originales cada componente principal:

\[\frac{\lambda_1}{\sum_{i=1}^{n} \lambda_i } \]

porc_explic=100*eigen1$values/sum(eigen1$values)
porc_explic
## [1] 54.54545 27.27273 18.18182

Claramente este summary esta mal

summary(pca1)
## Importance of components:
##                          PC1   PC2       PC3
## Standard deviation     1.225 1.225 4.532e-17
## Proportion of Variance 0.500 0.500 0.000e+00
## Cumulative Proportion  0.500 1.000 1.000e+00

c) Hallar los loadings de la primer componente principal.

Los loadings del primer componente principal son justamente el primer autovector que corresponde al autovalor mas grande, en este caso 6

Esto permite graficar los componentes, lo que hace es agregar un eje x de 1 a 3 (tiene dimensión 3)y luego crea barras para cada coeficiente usando geom_segment

df=data.frame(orden=as.factor(1:3),loadings=eigen1$vectors[,1],y=rep(0,3))
df
##   orden   loadings y
## 1     1 -0.4082483 0
## 2     2 -0.4082483 0
## 3     3 -0.8164966 0
ggplot(df )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,fill=orden,
                             color=orden,size=1.5))+ theme(legend.position = "none")

d) Hallar los scores de las primeras dos componentes principales correspondientes a la observación X=(2,2,1).

ACá la resolución de la profesora dice que no hay información para hacerlo porque faltan las medias, sin embargo el enunciado indica que tienen media 0. Entiendo que con ese dato si se puede hacer. Dejo la resolución para otro momento.

Ejercicio 2.

Considerando los datos de la base chalets.xls, se pide: ### SETUP

chalets <- read_excel("C:/Users/marco/Dropbox/Austral/Analisis_inteligente_de_datos/Clase_03/chalets.xls")

names(chalets) <- c("promotoras","hipoteca","precio_medio","sup_media_cocina")

a) Graficar el boxplot de cada una de las variables. Indicar, si se observa, la presencia de valores atípicos.

data_mod <- melt(chalets, id.vars='promotoras', 
                  measure.vars=c("hipoteca","precio_medio","sup_media_cocina"))

ggplot(data_mod) +
geom_boxplot(aes(x=promotoras, y=value, fill=variable,alpha=0.2)) +
  theme_minimal()

boxplot(chalets$hipoteca, plot=FALSE)$out
## [1] 37.3
boxplot(chalets$precio_medio, plot=FALSE)$out
## numeric(0)
boxplot(chalets$sup_media_cocina, plot=FALSE)$out
## numeric(0)

Vemos que el únio outlier es el de la hipoteca, las otras dos variables no tienen outliers

b) Graficar los diagramas de dispersión de las variables de a pares. Estimar la presencia de correlación entre variables a partir de estos gráficos, indicando si le parece fuerte y el signo de las mismas.

ggpairs(chalets,                 # Data frame
        columns = 2:4,        # Columns
        aes(alpha = 0.5))     # Transparency

c) Calcular el vector de medias y la matriz de varianzas y covarianzas muestral.

chalets2 <- as_tibble(t(colMeans(chalets[2:4])))  
chalets2 %>%   kbl(caption = "Vector de medias") %>% kable_classic(full_width =F)
Vector de medias
hipoteca precio_medio sup_media_cocina
19.05 1.57 9.73 
covb <- cov(chalets[2:4])
 covb %>%   kbl(caption = "Matriz de covarianzas") %>% kable_classic(full_width =F)
Matriz de covarianzas
hipoteca precio_medio sup_media_cocina
hipoteca 63.29833 5.7450000 33.863889
precio_medio 5.74500 0.9934444 4.053222
sup_media_cocina 33.86389 4.0532222 20.849000

d) Hallar la matriz de correlación muestral. Verificar las estimaciones realizadas visualmente.

``

corb <- cor(chalets[2:4])
 corb %>%   kbl(caption = "Matriz de correlaciones") %>% kable_classic(full_width =F)
Matriz de correlaciones
hipoteca precio_medio sup_media_cocina
hipoteca 1.0000000 0.7244728 0.9321763
precio_medio 0.7244728 1.0000000 0.8906068
sup_media_cocina 0.9321763 0.8906068 1.0000000

e) A partir de estas observaciones, le parece razonable pensar en un análisis de componentes principales para reducir la dimensión del problema?.

Si porque las variables esta altamente correlacionadas

f) Hallar la primera componente principal y graficar sus coeficientes mediante barras verticales.

pca1 <- prcomp(chalets[2:4], center=TRUE,scale=TRUE)
pca1$rotation  %>%   kbl(caption = "COnmponente Principal") %>% kable_classic(full_width =F)
COnmponente Principal
PC1 PC2 PC3
hipoteca 0.5687284 0.6607949 -0.4897940
precio_medio 0.5583782 -0.7474108 -0.3599874
sup_media_cocina 0.6039552 0.0687553 0.7940471 
vcont = chalets[2:4]
df1=data.frame(orden=as.factor(colnames(vcont)),loadings=pca1$rotation[,1],y=rep(0,3))
df1
##                             orden  loadings y
## hipoteca                 hipoteca 0.5687284 0
## precio_medio         precio_medio 0.5583782 0
## sup_media_cocina sup_media_cocina 0.6039552 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

g) Indicar qué porcentaje de la variabilidad total logra explicar esta componente. Explicar si se trata de una componente de tamaño o de forma. Es posible ordenar las promotoras en función de esta componente?. Si la respuesta es afirmativa, cual es la mayor y cual la menor; si es negativa, explicar por qué no es posible ordenarlos

summary(pca1)
## Importance of components:
##                           PC1     PC2     PC3
## Standard deviation     1.6435 0.52684 0.14574
## Proportion of Variance 0.9004 0.09252 0.00708
## Cumulative Proportion  0.9004 0.99292 1.00000

Explica el 90% de la varianza por si sóla. Es una componente de tamaño porque todas las barras tienen el mismo signo.

scores=cbind(as.data.frame(summary(pca1)$x),Promotora=1:10)

scores%>%arrange(-PC1)
##            PC1          PC2          PC3 Promotora
## 1   3.04473023  0.794448425 -0.076087522         8
## 2   2.31667376 -0.727038298  0.002634479        10
## 3   0.88274681 -0.674929544  0.178597294         7
## 4   0.02157016 -0.195921355 -0.200783631         3
## 5  -0.46084036  0.601313047 -0.095959157         5
## 6  -0.47859083  0.568280112  0.223793952         4
## 7  -0.85599312 -0.040926259  0.067138478         6
## 8  -1.02308426  0.072814352  0.129217853         2
## 9  -1.11892537 -0.390916578 -0.171447361         9
## 10 -2.32828702 -0.007123901 -0.057104384         1

Ejercicio 3.

Dado el siguiente conjunto de datos:

\[ X = \begin{equation} \begin{pmatrix} 3 & 6 \\ 5 & 6 \\ 10 & 12 \end{pmatrix} \end{equation} \]

matEj3 <- matrix(data = c(3, 6, 5, 6,10,12), nrow = 3, ncol = 2)

a) Calcule la matriz de covarianza, los autovalores y autovectores.

covEj3 <- cov(matEj3)
eigen3 <- eigen(covEj3)
aVal3 <- eigen3$values
aVec3 <- eigen3$vectors

covEj3 %>%   kbl(caption = "Matriz de varianzas y covarianzas") %>% kable_classic(full_width =F)
Matriz de varianzas y covarianzas
2.333333 3.666667
3.666667 9.333333 
aVal3 %>%   kbl(caption = "Autovalores") %>% kable_classic(full_width =F)
Autovalores
x
10.9023021
0.7643646 
aVec3 %>%   kbl(caption = "Autovectores") %>% kable_classic(full_width =F)
Autovectores
0.3933982 -0.9193682
0.9193682 0.3933982

b) Las componentes principales y su contribución porcentual a la varianza total.

pca3 <- prcomp(matEj3, center=FALSE,scale=FALSE)
pca3$rotation  %>%   kbl(caption = "Conmponente Principal") %>% kable_classic(full_width =F)
Conmponente Principal
PC1 PC2
-0.4441161 -0.8959693
-0.8959693 0.4441161 
summary(pca3)$importance[2:3,] %>%   kbl(caption = "Proporción explicada") %>% kable_classic(full_width =F)
Proporción explicada
PC1 PC2
Proportion of Variance 0.99544 0.00456
Cumulative Proportion 0.99544 1.00000

c) Grafique los datos en \(\mathbb{R}^{2x2}\) en la base original y en la base de los dos primeros ejes.

No entiendo que pide

d) Repita los cálculos con los datos estandarizados. Interprete los resultados obtenidos

corEj3 <- cor(matEj3)
eigen3d <- eigen(corEj3)
aVal3d <- eigen3d$values
aVec3d <- eigen3d$vectors
pca3d <- prcomp(matEj3, center=TRUE,scale=TRUE)

No me queda claro que se espera de la interpretación. En este caso no parece haber pesos de variables que distorsionen la representacion de la variabilidad por la escala. En un caso u otro las componentes representan variabilidad similar.

e) Verifique que los dos primeros autovectores son ortogonales entre sí. Represente gráficamente estos dos vectores en un gráfico bidimensional y trace rectas desde el origen hasta la ubicación de cada uno de los vectores en el gráfico.

El producto escalar entre ellos tiene que ser 0

aVec3d[,1] %*% aVec3d[,2]
##      [,1]
## [1,]    0

No entiendo que pide para graficar

Ejercicio 4.

Sea S la matriz de varianzas y covarianzas poblacionales:

mat4 <- matrix(data = c(3, 1, 1, 1,4,0,1,0,2), nrow = 3, ncol = 3)
mat4
##      [,1] [,2] [,3]
## [1,]    3    1    1
## [2,]    1    4    0
## [3,]    1    0    2

Correspondiente al vector aleatorio X = (X1,X2,X3)′ donde:

X1 puntuación media obtenida en las asignaturas de econometría X2 puntuación media obtenida en las asignaturas de derecho X3 puntuación media obtenida en asignaturas libres

Los datos corresponden a un conjunto d alumnos de la carrera de economía.

a) Calcule los autovalores de la matriz S.

eigen4 <- eigen(mat4)
eigen4
## eigen() decomposition
## $values
## [1] 4.732051 3.000000 1.267949
## 
## $vectors
##            [,1]       [,2]       [,3]
## [1,] -0.5773503  0.5773503  0.5773503
## [2,] -0.7886751 -0.5773503 -0.2113249
## [3,] -0.2113249  0.5773503 -0.7886751

b) Interprete la segunda componente principal sabiendo que el autovector correspondiente: e1 = (0, 5744;−0, 5744; 0, 5744).

Esta es una componente de forma ya que no tiene toda el mismo signo. Pondera de manera similar a las asignaturas libres y la economía y penaliza el saber derecho. Esta componente es fuerte en alumnos con perfiles poco orientados a conocer los aspectos legales.

c) Como se debería interpretar si un estudiante tuviera una puntuación en la componente principal muy inferior a la de sus compañeros?.

Que debería haber estudiado derecho y no economía

d) ¿Cuántas componentes principales serán necesarias para explicar al menos el 80% de la variabilidad total del conjunto?

pca4 <- prcomp(mat4,scale=FALSE)
summary(pca4)
## Importance of components:
##                           PC1    PC2      PC3
## Standard deviation     2.3154 1.1427 1.37e-16
## Proportion of Variance 0.8041 0.1959 0.00e+00
## Cumulative Proportion  0.8041 1.0000 1.00e+00

Alcanza con las dos primeras componentes Pero me da diferente que en la resolución ¿por qué?

autovalores <- eigen4$values
contribucion<-autovalores/sum(autovalores)
rbind(contribucion,cumsum(contribucion))
##                   [,1]      [,2]      [,3]
## contribucion 0.5257834 0.3333333 0.1408832
##              0.5257834 0.8591168 1.0000000

Ejercicio 5.

El siguiente conjunto de datos de la tabla 1 se refiere a 20 observaciones de suelo, donde se midió:

x1: contenido de arena, x2: contenido de cieno, x3: contenido de arcilla, x4: contenido de materia orgánica x5: acidez, según PH.

Compare los resultados del Análisis en Componentes Principales para la matriz de covarianza y para la matriz de correlación.

suelos <- read_excel("suelos.xlsx")
suelos <- suelos  %>%rename(MatOrg="Materia orgánica")
suelos %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Propiedades del suelo") %>% autofit()

a) Los porcentajes de variabilidad que logran explicar cada una de las componentes son los mismos?.

pcaSuelos <- prcomp(suelos[,-1],scale=FALSE)
summary(pcaSuelos)
## Importance of components:
##                            PC1    PC2     PC3     PC4       PC5
## Standard deviation     14.9613 2.8667 0.68736 0.50838 3.368e-15
## Proportion of Variance  0.9616 0.0353 0.00203 0.00111 0.000e+00
## Cumulative Proportion   0.9616 0.9969 0.99889 1.00000 1.000e+00
pcaSuelosEsc <- prcomp(suelos[,-1],scale=TRUE)
summary(pcaSuelosEsc)
## Importance of components:
##                           PC1    PC2    PC3     PC4       PC5
## Standard deviation     1.6311 1.0713 0.9810 0.47891 4.107e-16
## Proportion of Variance 0.5321 0.2295 0.1925 0.04587 0.000e+00
## Cumulative Proportion  0.5321 0.7616 0.9541 1.00000 1.000e+00
#fviz_screeplot(pcaSuelos, addlabels = TRUE, ylim = c(0,97))
#fviz_pca_ind(pcaSuelos, geom.ind = "point", 
#             col.ind = "#FC4E07", 
#             axes = c(1, 2), 
#             pointsize = 1.5) 
fviz_contrib(pcaSuelos, choice = "var", axes = 1, top = 10)

#biplot(pcaSuelos, scale = 0, cex = 0.5, col = c("dodgerblue3", "deeppink3"))
fviz_contrib(pcaSuelosEsc, choice = "var", axes = 1, top = 10)

Si cambian los porcentajes que explica cada variable. La razón detás de esto es la diferencia de escala entre las variables originales, en el gráfico de contribuciones de cada variable original, podemos ver como pesa la variable arena en el caso de la matriz de covarianza. Cuando escalamos y centramos con la matriz de correlación esto cambia.

b) Cambia el orden de las componentes?

pcaSuelos$rotation
##                   PC1          PC2         PC3          PC4           PC5
## Arena   -0.7849449544 -0.223100113  0.02718830 -0.003901445 -5.773503e-01
## Cieno    0.5870812022 -0.560792619  0.08607941 -0.010212850 -5.773503e-01
## Arcilla  0.1978637522  0.783892732 -0.11326771  0.014114294 -5.773503e-01
## MatOrg   0.0068110819 -0.145758773 -0.97995170  0.135656360  7.806256e-18
## Acidez   0.0007907581 -0.002131353 -0.13680723 -0.990595081 -5.681219e-17
pcaSuelosEsc$rotation
##                 PC1         PC2         PC3         PC4          PC5
## Arena   -0.61054577  0.02208448 -0.03234998 -0.17078162 7.723557e-01
## Cieno    0.58144760 -0.14575518  0.10497554  0.53458750 5.864047e-01
## Arcilla  0.53495413  0.28024523 -0.14981128 -0.74380575 2.441235e-01
## MatOrg   0.04931997 -0.83603278  0.40937786 -0.36197582 1.249001e-16
## Acidez   0.02332740 -0.44808426 -0.89325246  0.02786334 2.775558e-17

mmm… son distintos….cambian, no se si el orden…

c) Cambian los loadings de las componentes?

Otra función que calcula PCA…. en este caso vemos los loadings

suelos.pca.corr <- princomp(suelos[,-1], cor = FALSE,scores=TRUE)
suelos.pca.corr$loadings
## 
## Loadings:
##         Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Arena    0.785  0.223                0.577
## Cieno   -0.587  0.561                0.577
## Arcilla -0.198 -0.784 -0.113         0.577
## MatOrg          0.146 -0.980 -0.136       
## Acidez                -0.137  0.991       
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings       1.0    1.0    1.0    1.0    1.0
## Proportion Var    0.2    0.2    0.2    0.2    0.2
## Cumulative Var    0.2    0.4    0.6    0.8    1.0
suelos.pca.corr <- princomp(suelos[,-1], cor = TRUE,scores=TRUE)
suelos.pca.corr$loadings
## 
## Loadings:
##         Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Arena    0.611                0.171  0.772
## Cieno   -0.581 -0.146 -0.105 -0.535  0.586
## Arcilla -0.535  0.280  0.150  0.744  0.244
## MatOrg         -0.836 -0.409  0.362       
## Acidez         -0.448  0.893              
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings       1.0    1.0    1.0    1.0    1.0
## Proportion Var    0.2    0.2    0.2    0.2    0.2
## Cumulative Var    0.2    0.4    0.6    0.8    1.0

Los loadings cambian porque el peso relativo de cada variable cambia en los componentes y los autovectores resultantes en cada caso son diferentes.

d) Cuál de los dos análisis le parece más adecuado y por qué?.

Siempre es más robusto trabajar con la matríz de correlaciones para evitar los problemas de escala.

Ejercicio 6. La tabla gorriones.xls contiene datos de 49 aves, 21 de los cuales sobrevivieron a una tormenta. Se pide:

gorriones <- read_excel("gorriones.xlsx")
head(gorriones) %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Primeros valores de gorriones.xlsx") %>% autofit()

a) Estandarice las variables y calcule la matriz de covarianzas para las variables estandarizadas.

as.data.frame(cov(scale(gorriones[,-1]))) %>% add_rownames(" ") %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Matriz de varianza y covarianza de matriz estandarizada") %>% autofit()
## Warning: `add_rownames()` was deprecated in dplyr 1.0.0.
## ℹ Please use `tibble::rownames_to_column()` instead.

b) Verifique que ésta es la matriz de correlación de las variables originales.

corr <- cor(gorriones[,-1])
as.data.frame(corr)  %>% add_rownames(" ") %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Matriz de correlación") %>% autofit()

Verificado ### c) Le parece adecuado en este caso un análisis de componentes principales. ¿Qué indica el autovalor para una componente principal? De la tabla de correlaciones se puede ver que hay varias variables con fuerte correlación, en este caso aplica el análisis de componentes principales

ggcorrplot(as.data.frame(corr))

d) ¿Cuántas componentes son necesarias para explicar el 80% de la varianza total? Realice el grafico de sedimentación, fundamente su respuesta con este gráfico.

pcaGorr <- prcomp(gorriones[,-1], scale=TRUE)
sumPcaGorr <- summary(pcaGorr)
as.data.frame(sumPcaGorr$importance) %>% add_rownames(" ") %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Importancia de las componentes principales") %>% autofit()
fviz_screeplot(pcaGorr, addlabels = TRUE, ylim = c(0,97))

En el gráfico se ve que a partir del tercer componente la variabilidad que explican los siguientes es muy poca, dependiendo el corte que se desee hacer podemos quedarnos con 3 o cuatro componentes. Si usamoes el criterio de Kaiser tomariamos todos los componentes cuyos respectivos autovalores \(\lambda_m \geq 1\) y agregarle \(\lambda_{m+1} \leq 1\) (o \(\lambda_m \geq 0.7\) en la forma ajustada) si trabajamos con la matriz de correlación

data.frame(Autovalores=eigen(corr)$values) %>% flextable() %>% align_text_col(align = "right") %>% autofit()

Sin embargo el test de Kaiser nos dice que hay que tomar sólo las dos primeras si tomamos el criterio de \(\lambda_m \geq 0.7\) y no el de \(\lambda_m \geq 1\) + \(\lambda_{m+1} \leq 1\). En realidad son las que más explican, pero la tercera tiene algo que decir y según el criterio del gráfico debería considerarse

e) ¿Cuál es la expresión de la primer componente principal?

El primer componente principal es el que usa el primer autovector como coordenadas y las variables medidas normalizadas.

as.data.frame(pcaGorr$rotation) %>% add_rownames(" ") %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Componentes principales") %>% autofit()

entonces la primer componente sería:

\[Y_1=0,46*Z_1+0,46*Z_2+0.45*Z_3+0.447*Z_4+0,40*Z_5+-0,03*Z_7\] ### f) ¿Cómo queda expresada la primer componente principal? (en función del autovector correspondiente y de las variables).

¿Cuál sería la diferencia con el punto anterior?

g) Encuentre las coordenadas del pájaro 11 en las nuevas componentes.

Es el producto vectorial del vector de coordenadas y los valores para el pajaro 11 pero estandarizadas.

gorrScaled=scale(gorriones[,2:7])

data.frame(
  PC1= as.vector(pcaGorr$rotation[,1]) %*% as.vector(t(gorrScaled[11,])),
  PC2= as.vector(pcaGorr$rotation[,2]) %*% as.vector(t(gorrScaled[11,])),
  PC3= as.vector(pcaGorr$rotation[,3]) %*% as.vector(t(gorrScaled[11,]))
  )
##         PC1      PC2      PC3
## 1 0.3611794 1.193108 1.151008

Acá lo hacemos para todas las componentes

as.vector(gorrScaled[11,]) %*% as.matrix(pcaGorr$rotation) 
##            PC1      PC2      PC3        PC4        PC5       PC6
## [1,] 0.3611794 1.193108 1.151008 0.01290467 -0.0328486 0.1137724

h) Represente gráficamente en el plano. (Eje 1 vs 2, 1 vs 3, 2 vs 3). Interprete los tres primeros ejes.

fviz_pca_biplot(pcaGorr, axes=c(1,2), repel = TRUE)

fviz_pca_biplot(pcaGorr, axes=c(1,3), repel = TRUE)

fviz_pca_biplot(pcaGorr, axes=c(2,3), repel = TRUE)

Agrego la caracterización de las dimensiones para tener como referencia:

df1=data.frame(orden=as.factor(colnames(gorriones[2:7])),loadings=pcaGorr$rotation[,1],y=rep(0,6))
df1
##                 orden   loadings y
## largototal largototal  0.4557344 0
## extension   extension  0.4602511 0
## cabeza         cabeza  0.4476819 0
## humero         humero  0.4716751 0
## esternon     esternon  0.3960364 0
## sobrevida   sobrevida -0.0270888 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

df1=data.frame(orden=as.factor(colnames(gorriones[2:7])),loadings=pcaGorr$rotation[,2],y=rep(0,6))
df1
##                 orden    loadings y
## largototal largototal -0.14676423 0
## extension   extension -0.02017345 0
## cabeza         cabeza  0.03012167 0
## humero         humero  0.16357298 0
## esternon     esternon  0.03011804 0
## sobrevida   sobrevida  0.97441406 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

df1=data.frame(orden=as.factor(colnames(gorriones[2:7])),loadings=pcaGorr$rotation[,3],y=rep(0,6))
df1
##                 orden     loadings y
## largototal largototal  0.014596454 0
## extension   extension -0.306278675 0
## cabeza         cabeza -0.312035697 0
## humero         humero -0.161512413 0
## esternon     esternon  0.884590279 0
## sobrevida   sobrevida  0.005274471 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

### i) Realice un gráfico donde se observen los gorriones en los nuevos ejes 1 y 2, y resalte con distinto color el grupo de los que sobrevivieron.

fviz_pca_ind(pcaGorr, axes=c(1,2), repel = TRUE, habillage = gorriones$sobrevida, addEllipses=TRUE, ellipse.level=0.95)

j) Utilice el Análisis en Componentes Principales como método para encontrar outliers.

NO recuerdo esto

Ejercicio 7. Con el objetivo de oObtener índices útiles para la gestión hospitalaria basados en técnicas estadísticas multivariantes descriptivas se se recogió información del Hospital de Algeciras correspondiente a los ingresos hospitalarios del periodo 2007-

Se estudiaron las variables habitualmente monitorizadas por el Servicio Andaluz de Salud, del Sistema Nacional de Salud Español:

  • NI: número de ingresos
  • MO: mortalidad
  • RE: número de reingresos
  • NE: número de consultas externas
  • ICM: Indice cardíaco máximo
  • ES: número de estancias
  • IF: Variable que agregaron para el ejercicio que no se que sería

Las variables se midieron en un total de 22486 ingresos.

En la siguiente tabla se aprecia la Distribución de los valores obtenidos en las variables listadas por los servicios del hospital de Algeciras( Andalucía,Espana):

hospitales <- read_excel("hospitales.xlsx")
hospitales %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Por servicio del Hospital") %>% autofit()

La idea central del ACP es conseguir la simplificación de un conjunto de datos, generalmente cuantitativos, procedentes de un conjunto de variables interrelacionadas. Este objetivo se alcanza obteniendo, a partir de combinaciones lineales de las variables originalmente medidas, un nuevo conjunto de igual número de variables, no correlacionadas, llamadas componentes principales (CP) en las cuales permanece la variabilidad presente en los datos originales, y que al ordenarlas decrecientemente por su varianza, nos permiten explicar el fenómeno de estudio con las primeras CP.

Verificar que las primeras dos componentes principales son:

\[Y_1 = 0,5380NI+0,5126ES+0,4081IF+0,2635MO−0,1561NE−0,2535RE− 0,3511ICM\]

\[Y_2 = 0,5524MO+0,4952RE+0,4696ICM+0,3756ES+0,2867NE+0,05778IF− 0,04908NI\]

a) Grafique las cargas y explicar la interpretación de las componentes principales.

pcaHosp <- prcomp(hospitales[2:8],scale=TRUE)

as.data.frame(pcaHosp$rotation)  %>% add_rownames(" ") %>% flextable() %>% align_text_col(align = "right") %>% set_caption(caption = "Componentes Principales") %>% autofit()

Esto verifica los componentes indicados.

Acá me entra la duda, en teoría los autovectores calculados sobre la matríz de correlación, deberían aplicarse sobre las variables estandarizadas, sin embargo el enunciado no parece decir eso, perop las coordenadas coinciden….

sumHosp <- summary(pcaHosp)

as.data.frame(sumHosp$importance[2:3,])  %>% add_rownames(" ") %>% flextable()%>% align_text_col(align = "right") %>% set_caption(caption = "Contribución a la  explicación de la varianza") %>% autofit()

Estos son los autovalores correspondientes wl paquete factorMineR con la función PCA los devuelve como parametro, hay que explorar más esa función y ese paquete

eigen(cor(hospitales[2:8]))$values
## [1] 2.55803188 1.82896471 1.24330111 0.64842333 0.51187926 0.16302704 0.04637267

Si seguimos el criterio de Kaiser deberiamos quedarnos con los 4 primeros componentes.

Graficamos entonces estas 4 componentes:

df1=data.frame(orden=as.factor(colnames(hospitales[2:8])),loadings=pcaHosp$rotation[,1],y=rep(0,7))
df1
##     orden   loadings y
## NI     NI -0.5378660 0
## MO     MO -0.2635670 0
## RE     RE  0.2535153 0
## NE     NE  0.1564391 0
## ICM   ICM  0.3510125 0
## ES     ES -0.5125446 0
## IF     IF -0.4081393 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

df1=data.frame(orden=as.factor(colnames(hospitales[2:8])),loadings=pcaHosp$rotation[,2],y=rep(0,7))
df1
##     orden    loadings y
## NI     NI  0.04913727 0
## MO     MO -0.55250078 0
## RE     RE -0.49523090 0
## NE     NE -0.28625361 0
## ICM   ICM -0.46987734 0
## ES     ES -0.37550358 0
## IF     IF -0.05784453 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

df1=data.frame(orden=as.factor(colnames(hospitales[2:8])),loadings=pcaHosp$rotation[,3],y=rep(0,7))
df1
##     orden   loadings y
## NI     NI  0.3242016 0
## MO     MO -0.1945829 0
## RE     RE  0.2053858 0
## NE     NE  0.6486612 0
## ICM   ICM -0.4073489 0
## ES     ES  0.1432844 0
## IF     IF -0.4556562 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

df1=data.frame(orden=as.factor(colnames(hospitales[2:8])),loadings=pcaHosp$rotation[,4],y=rep(0,7))
df1
##     orden    loadings y
## NI     NI -0.37373922 0
## MO     MO  0.35580106 0
## RE     RE -0.65347352 0
## NE     NE  0.53974190 0
## ICM   ICM -0.05946509 0
## ES     ES -0.06146428 0
## IF     IF  0.08978696 0
ggplot(df1 )+geom_segment(aes(x=orden,xend=orden,y=y,yend=loadings,color=orden,
                             size=1.5))+ theme(legend.position = "none")+xlab("")

### b) Qué porcentaje de variabilidad logra captar cada una de ellas?. Grafique el scree plot.

fviz_screeplot(pcaHosp, addlabels = TRUE, ylim = c(0,97))

Como se menciono en a) se ve que hay que tomar las 4 primeras componentes. ### c) Le parece adecuado considerar dos componentes principales?.

Contestado en a), no me parece que por criterio de kaiser debemos quedarnos con 4. ### d) Hallar la correlación entre las nuevas variables y las originales.

M = cor(hospitales[2:8])
corrplot(M)

Calculo las nuevas variables con las coordenadas de los autovectores aplicadas a las variables medidas estandarizadas

scores <- as.matrix(scale(hospitales[2:8])) %*% as.matrix(pcaHosp$rotation)
M = cor(scores)
corrplot(M)

Y comprobamos que todos los componentes son ortogonales, es decir no tienen correlación entre ellos.

e) Ordenar los servicios en función de su puntuación en cada una de las dos primeras componentes principales. Indicar cuales son los servicios más demandados y los más complejos.

scoresDF <- as.data.frame(scores) 
rownames(scoresDF) <- t(hospitales[,1])
scoresDF %>% arrange(desc(PC1))  %>% add_rownames(" ") %>% flextable()%>% align_text_col(align = "right") %>% set_caption(caption = "Contribución a la  explicación de la varianza") %>% autofit()

No sabría como interpretar complejidad en base a los datos del ejercicio, pero la primera componente favorece reingresos e indice cardiaco medio, podriamos decir que es esa. O también la componente 2 con signo cambiado donde penaliza los ingresos…. Quedemosnos con esa….yo que se…

En cuanto a los mas demandados probablemente sean los que tengan elevada la \(3^{ra}\) componente.

f) Representar un biplot y buscar servicios similares, asociaciones entre las variables. Verificar en este grafico la representación de las variables originales en las componentes.

fviz_pca_biplot(pcaHosp)

  • La especialidad psiquiatrica es la que mas puntúa en la primera componente. Más reingresos e ICM. En el lado opuesto esta medicina interna

*12 y 13 son similares y pero cercanos al promedio. 2 y 10 son similares y 7 y 3 también.