Resolución Examen M4
Borja Macías Afonso
3 de junio de 2021
library(readxl)
library(foreign)
library(vegan)
library(smacof)
library(FactoMineR)
library(factoextra)
library(dendextend)
library(MASS)
library(tidyverse)
library(lubridate)
library(patchwork)
library(formattable)
library(gridExtra)
library(corrplot)
library(ggpubr)
library(psych)EJERCICIO 1
datos_covid<- read_xlsx("datos_covid.xlsx")
str(datos_covid)## tibble [9,234 x 8] (S3: tbl_df/tbl/data.frame)
## $ ccaa_iso : chr [1:9234] "AN" "AR" "AS" "CB" ...
## $ fecha : POSIXct[1:9234], format: "2020-01-01" "2020-01-01" ...
## $ num_casos : num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...
## $ num_casos_prueba_pcr : num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...
## $ num_casos_prueba_test_ac : num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...
## $ num_casos_prueba_ag : num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...
## $ num_casos_prueba_elisa : num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...
## $ num_casos_prueba_desconocida: num [1:9234] 0 0 0 0 0 0 0 0 0 0 ...En primer lugar, realizamos las modificaciones necesarias para obtener una tabla de frecuencias (número de casos por meses y CCAA) que nos permita calcular las distancias entre CCAA. Para ello, creamos una nueva variable que señala sólamente el mes y el año del registro (al día le pone siempre el valor 1). Posteriormente, agrupamos los registros por CCAA y meses sumando en número de casos y pivotamos la tabla para que cada mes sea una variable (columna). Finalmente, pasamos la columna del nombre de las comunidades a nombre de filas.
datos_covid$anio_mes<-make_date(year(datos_covid$fecha),month(datos_covid$fecha)) #Nueva variable indicando sólo año y mes
datos_covid2<- group_by(datos_covid,ccaa_iso,anio_mes) %>% summarise(casos=sum(num_casos)) %>% pivot_wider(names_from = anio_mes,values_from=casos) %>% column_to_rownames("ccaa_iso")
formattable(head(datos_covid2))| 2020-01-01 | 2020-02-01 | 2020-03-01 | 2020-04-01 | 2020-05-01 | 2020-06-01 | 2020-07-01 | 2020-08-01 | 2020-09-01 | 2020-10-01 | 2020-11-01 | 2020-12-01 | 2021-01-01 | 2021-02-01 | 2021-03-01 | 2021-04-01 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AN | 0 | 73 | 10038 | 2778 | 326 | 446 | 2649 | 20426 | 34741 | 90433 | 77072 | 31903 | 160468 | 50763 | 28275 | 43924 |
| AR | 0 | 23 | 3484 | 2018 | 432 | 546 | 8671 | 11989 | 10096 | 22884 | 13737 | 7128 | 18647 | 7476 | 4305 | 7689 |
| AS | 0 | 5 | 1593 | 752 | 76 | 11 | 102 | 941 | 2290 | 7469 | 10946 | 3355 | 10955 | 5470 | 3737 | 3294 |
| CB | 0 | 11 | 1509 | 714 | 85 | 48 | 161 | 2045 | 2415 | 3724 | 4926 | 2605 | 4642 | 1943 | 1425 | 2737 |
| CE | 0 | 0 | 88 | 23 | 42 | 10 | 3 | 158 | 388 | 1290 | 847 | 274 | 931 | 582 | 587 | 536 |
| CL | 0 | 60 | 12280 | 7127 | 989 | 554 | 946 | 10089 | 20680 | 40133 | 31512 | 9635 | 58286 | 14561 | 6928 | 10080 |
Para calcular las distancias entre CCAA utilizaremos frecuencias acumuladas en lugar de las frecuencias puntuales de cada mes, ya que dan una mejor idea de la evolución de los casos.
acumulado<- datos_covid2 %>% apply(MARGIN = 1,cumsum)
acumulado<-as.data.frame(t(acumulado))
formattable(slice_sample(acumulado)) # Ejemplo de comprobación| 2020-01-01 | 2020-02-01 | 2020-03-01 | 2020-04-01 | 2020-05-01 | 2020-06-01 | 2020-07-01 | 2020-08-01 | 2020-09-01 | 2020-10-01 | 2020-11-01 | 2020-12-01 | 2021-01-01 | 2021-02-01 | 2021-03-01 | 2021-04-01 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CT | 2 | 249 | 34145 | 54743 | 59738 | 63344 | 85366 | 115718 | 147338 | 248860 | 310044 | 366054 | 464818 | 507044 | 542996 | 587614 |
Para seleccionar el tipo de distancia a calcular hay que tener en cuenta dos aspectos:
- El tipo de datos. En este caso estamos ante series temporales de frecuencias absolutas.
- El objetivo del análisis. Lo que queremos evaluar es si existen diferencias en la “evolución” de los casos por CCAA. Para ello, es necesario utilizar medidas que tengan en cuenta la masa o escala de cada CCAA, porque de lo contrario lo único que haríamos con el escalamiento es ordenar las comunidades de mayor a menor número de casos totales (estando esto condicionado por la cantidad de población en cada una), y este no es el objetivo. Lo que queremos es saber si existen diferencias en la distribución del total de casos de cada CCAA a lo largo de los 16 meses de estudio.
En base a estos criterios calcularemos la distancia chi-cuadrado y el coeficiente de correlación de Pearson (transformado a distancia).
# Distancia chi-cuadrado
dist_chisq<-vegdist(acumulado,method = "chisq")
# Distancia de correlación
correl<-cor(t(acumulado))
dist_correl<-sqrt(2*(1-correl))Para seleccionar una de las dos medidas para realizar el escalamiento escogeremos la que tenga menor stress. Probaremos tanto con escalamiento métrico como no métrico:
tipo<-c("ratio","interval","ordinal")
distancias<-list(dist_chisq,dist_correl)
nombre<-c("chi-cuadrado","correlación")
for(i in 1:length(distancias)){
for(j in 1:length(tipo)){
d<-distancias[[i]]
escalamiento<- mds(d,type = tipo[j])
print(nombre[i])
print(tipo[j])
print(escalamiento$stress)
}}## [1] "chi-cuadrado"
## [1] "ratio"
## [1] 0.09777422
## [1] "chi-cuadrado"
## [1] "interval"
## [1] 0.07605097
## [1] "chi-cuadrado"
## [1] "ordinal"
## [1] 0.06199794
## [1] "correlación"
## [1] "ratio"
## [1] 0.1300671
## [1] "correlación"
## [1] "interval"
## [1] 0.1035518
## [1] "correlación"
## [1] "ordinal"
## [1] 0.08991438A la vista de los resultados optamos por un escalamiento no métrico con la distancia chi-cuadrado, que cuenta con un nivel bueno de stress de 6,2%.
esc_chisq<- mds(dist_chisq,type = "ordinal",ndim = 2)
ggplot()+ geom_point(aes(esc_chisq$conf[,1],esc_chisq$conf[,2]),col="red")+
geom_text(aes(esc_chisq$conf[,1],esc_chisq$conf[,2],label=rownames(datos_covid2),vjust=-0.5))+
labs(title="Escalamiento no métrico con chi-cuadrado",x="Dim 1",y="Dim 2")+
geom_vline(xintercept = 0,linetype="dashed")+
geom_hline(yintercept=0,linetype="dashed")+
theme_light()
Para comprender mejor lo que explica cada una de las dos dimensiones, graficaremos las series de las cuatro CCAA más extremas de cada cuadrante: La rioja, Valencia, Navarra y Melilla. Utilizaremos valores relativos para tener en cuenta la masa de cada CCAA.
# Obtención de valores relativos
totales<-data.frame(total=rowSums(datos_covid2)) %>% rownames_to_column(var = "ccaa_iso")
datos_prop<- data.frame(t(datos_covid2 /totales$total))
acumulado_prop<-as.data.frame(apply(datos_prop, 2,cumsum))
# Gráfico con datos puntuales
g11<- datos_prop %>%
rownames_to_column("anio_mes") %>%
select(anio_mes,RI,VC,NC,ML) %>%
pivot_longer(cols=c("RI","VC","NC","ML"),names_to = "ccaa",values_to="casos") %>%
ggplot()+
geom_line(mapping = aes(x=anio_mes,y=casos,group=ccaa,col=ccaa))+
labs(title="Puntuales")+
theme(axis.text.x = element_text(angle = 45,size = 7))
# Gráfico con datos acumulados
g12<- acumulado_prop %>%
rownames_to_column("anio_mes") %>%
select(anio_mes,RI,VC,NC,ML) %>%
pivot_longer(cols=c("RI","VC","NC","ML"),names_to = "ccaa",values_to="casos") %>%
ggplot()+
geom_line(mapping = aes(x=anio_mes,y=casos,group=ccaa,col=ccaa))+
labs(title="Acumulados")+
theme(axis.text.x = element_text(angle = 45,size = 7))
g11/g12
En el primer gráfico sobre el número de casos puntuales se pueden observar claramente las 3 olas de contagios que hemos vivido durante la pandemia. A la vista de estos dos gráficos, podemos concluir que la primera dimensión del escalamiento recoge la gravedad de la primera ola, puesto confronta comunidades muy afectadas por la primera ola situándolas a la izquierda (La Rioja, Navarra) con otras menos afectadas situándolas a la derecha (Comunidad Valenciana y Melilla). Por su parte, la segunda dimensión parece recoger la gravedad de la tercera ola, situando más arriba a las comunidades más afectadas por esta ola (Comunidad Valenciana y La rioja) y abajo las menos afectadas (Navarra y Melilla). Hay que recordar que hablamos de gravedad en términos relativos de casos (respecto al total de cada CCAA) y no en términos absolutos.
Podríamos añadir una dimensión más al escalamiento y probablemente reduciríamos el stress y recogería el efecto de la segunda ola, aunque la interpretación a la hora de buscar agrupaciones entre CCAA sería más complicada.
esc_3d<- mds(dist_chisq,type = "ordinal",ndim = 3)
esc_3d$stress## [1] 0.01738323# Gráficos de coordenadas para cada par de dimensiones
g13<-ggplot()+ geom_point(aes(esc_3d$conf[,1],esc_3d$conf[,2]),col="red")+
geom_text(aes(esc_3d$conf[,1],esc_3d$conf[,2],label=rownames(acumulado),vjust=-0.5),cex=2.5)+
labs(x="Dim 1",y="Dim 2")+
geom_vline(xintercept = 0,linetype="dashed")+
geom_hline(yintercept=0,linetype="dashed")+
theme_light()
g14<-ggplot()+ geom_point(aes(esc_3d$conf[,1],esc_3d$conf[,3]),col="red")+
geom_text(aes(esc_3d$conf[,1],esc_3d$conf[,3],label=rownames(acumulado),vjust=-0.5),cex=2.5)+
labs(x="Dim 1",y="Dim 3")+
geom_vline(xintercept = 0,linetype="dashed")+
geom_hline(yintercept=0,linetype="dashed")+
theme_light()
g15<-ggplot()+ geom_point(aes(esc_3d$conf[,2],esc_3d$conf[,3]),col="red")+
geom_text(aes(esc_3d$conf[,2],esc_3d$conf[,3],label=rownames(acumulado),vjust=-0.5),cex=2.5)+
labs(x="Dim 2",y="Dim 3",cex=1)+
geom_vline(xintercept = 0,linetype="dashed")+
geom_hline(yintercept=0,linetype="dashed")+
theme_light()
# Gráfico con datos puntuales
g16<-datos_prop %>%
rownames_to_column("anio_mes") %>%
select(anio_mes,IB,CN,AS,CE) %>%
pivot_longer(cols=c("IB","CN","AS","CE"),names_to = "ccaa",values_to="casos") %>%
ggplot()+
geom_line(mapping = aes(x=anio_mes,y=casos,group=ccaa,col=ccaa))+
labs(title="Puntuales")+
theme(axis.text.x = element_text(angle = 45,size = 7))
# Gráfico con datos acumulados
g17<-acumulado_prop %>%
rownames_to_column("anio_mes") %>%
select(anio_mes,IB,CN,AS,CE) %>%
pivot_longer(cols=c("IB","CN","AS","CE"),names_to = "ccaa",values_to="casos") %>%
ggplot()+
geom_line(mapping = aes(x=anio_mes,y=casos,group=ccaa,col=ccaa))+
labs(title="Acumulados")+
theme(axis.text.x = element_text(angle = 45,size = 7))
(g13|g14|g15)/(g16|g17)
El stress se reduce hasta un 1,7%, lo que indica una representatividad muy buena de la nueva matriz de distancias.
En cuanto a la interpretación de la tercera dimensión, parece separar las CCAA en función de la gravedad y momento en que se produce la segunda ola. Así, Canarias y Baleares se sitúan más arriba al ser las comunidades donde antes se vive la segunda ola. En contra, Asturias y Ceuta experimentan la segunda ola de forma más tardía, como se puede ver en el gráfico de número de casos puntuales por mes.
EJERCICIO 2
Realizamos en primer lugar un clúster jerarquico, para tener una primera idea de agrupación:
cluster<-hclust(esc_chisq$confdist,method = "ave")
plot(cluster,hang=-1)
Para seleccionar el número óptimo de clusters vamos a realizar distintas simulaciones con hasta 12 clusters y analizaremos la reducción del error intragrupos mediante un gráfico de sedimentación y también el tamaño de los clusters.
K_Max <- 12
Errores <- NULL
set.seed(12345)
# Análisis del número de clusters óptimo
for (i in 1:K_Max){
Errores[i] <- sum(kmeans(esc_chisq$conf, centers=i)$withinss)
if (i == 1){
Errores_max <- Errores[i]
}
Errores[i] <- Errores[i] / Errores_max
}
ggplot(mapping = aes(2:K_Max,Errores[-1]))+geom_point()+geom_line() # En base a la disminución del error intragrupos
for(i in 2:K_Max){
nclust<-cutree(cluster,k=i) #En base al número de observaciones en cada cluster
print(table(nclust))
}## nclust
## 1 2
## 9 10
## nclust
## 1 2 3
## 8 10 1
## nclust
## 1 2 3 4
## 6 10 2 1
## nclust
## 1 2 3 4 5
## 6 3 7 2 1
## nclust
## 1 2 3 4 5 6
## 6 3 4 2 3 1
## nclust
## 1 2 3 4 5 6 7
## 4 3 2 4 2 3 1
## nclust
## 1 2 3 4 5 6 7 8
## 4 3 2 3 2 3 1 1
## nclust
## 1 2 3 4 5 6 7 8 9
## 4 1 2 3 2 3 1 2 1
## nclust
## 1 2 3 4 5 6 7 8 9 10
## 4 1 2 3 2 2 1 2 1 1
## nclust
## 1 2 3 4 5 6 7 8 9 10 11
## 4 1 2 3 2 2 1 1 1 1 1
## nclust
## 1 2 3 4 5 6 7 8 9 10 11 12
## 2 1 2 3 2 2 2 1 1 1 1 1En base a la información anterior, realizaremos 6 clusters. Con dos clusters los grupos tienen una tamaño grande y homogeneo, sin embargo, el error intragrupos es demasiado alto (33% aprox.).Con 4 clusters ya conseguiamos un error intragrupos bueno por debajo del 20%, aunque hay dos clusters con solo 1 y 2 observaciones. Si pasamos de 4 a 6 clusters reducimos el error hasta un 0,1 aproximadamemte y el cluster de mayor tamaño(10) se desagrega en 3 de tamaño 4,3 y 3, por lo que los tamaños se hacen más homogeneos.
Así, el dendograma quedaría de la siguiente forma:
dend<-color_branches(cluster,k=6)
plot(dend)
Representando los grupos en un gráfico de coordenadas con las dos dimensiones obtenidas en el escalamiento y teniendo en cuenta las interpretaciones dadas en el ejercicio anterior podemos tener una idea más clara de la evolución de los diferentes grupos.
grupo<-cutree(cluster,k=6)
datos<-data.frame(esc_chisq$conf)
colnames(datos)<-c("Dim.1","Dim.2")
datos$grupo<-factor(grupo)
ggscatter(datos,x = "Dim.1", y = "Dim.2",
label = rownames(datos),
color = "grupo",
palette = "jco",
size = 1,
ellipse = TRUE,
ellipse.type = "convex",
repel = TRUE) +
geom_vline(xintercept = 0,linetype="dashed") +
geom_hline(yintercept=0,linetype="dashed")
Al gráfico anterior añadimos otro que ilustra la evolución de los diferentes grupos, obtenido a partir de las medias de casos por mes de cada grupo (nuevamente en términos relativos).
# Calculo de medias de cada grupo en cada mes
g21<-as.data.frame(t(datos_prop)) %>% rownames_to_column() %>% mutate(grupo=grupo) %>% group_by(grupo) %>% summarise(across(c(-rowname), mean)) %>% pivot_longer(cols=c(-grupo),names_to="anio_mes",values_to="casos") %>%
ggplot()+ # Gráfico de medias por grupo en cada mes
geom_line(mapping = aes(x=anio_mes,y=casos,group=grupo,col=as.factor(grupo)))+
labs(title="Puntuales")+
theme(axis.text.x = element_text(angle = 45,size = 7))
tema<-ttheme_default(
core = list(fg_params=list(cex = 0.8)),
colhead = list(fg_params=list(cex = 0.8)),
rowhead = list(fg_params=list(cex = 0.8))) # Para ajustar el tamaño del texto de la tabla
tabla<-tableGrob(datos_covid2 %>% rownames_to_column(var="ccaa") %>% mutate(grupo=grupo) %>% select(ccaa,grupo) %>% arrange(grupo) %>% t(),theme=tema) # Tabla de grupos
g21/tabla
A partir de los dos gráficos anteriores podemos describir la evolución de cada grupo de la siguiente forma:
- Grupo 1: se trata de comunidades con una incidencia muy baja en la primera ola, alta en la segunda y media durante la tercera.
- Grupo 2: incidencia media-alta durante la primera ola, alta en la segunda y muy baja en la tercera.
- Grupo 3: incidencia media-alta en la primera ola, media en la segunda y media-alta en la tercera.
- Grupo 4: incidencia baja en la primera ola, media en la segunda ola (más temprana que el resto de CCAA) y media en la tercera.
- Grupo 5: incidencia muy alta en la primera ola, media en la segunda y media en la tercera
- Grupo 6: incidencia baja en la primera ola, baja en la segunda y muy alta en la tercera.
EJERCICIO 3
En primer lugar, vamos a obtener la tabla de contingencia:
tipo_prueba<-datos_covid %>% group_by(ccaa_iso) %>% summarise(pcr=sum(num_casos_prueba_pcr),test_ac=sum(num_casos_prueba_test_ac),antigenos=sum(num_casos_prueba_ag),elisa=sum(num_casos_prueba_elisa),desconocido=sum(num_casos_prueba_desconocida)) %>% column_to_rownames(var="ccaa_iso")
formattable(tipo_prueba)| pcr | test_ac | antigenos | elisa | desconocido | |
|---|---|---|---|---|---|
| AN | 378442 | 323 | 175543 | 7 | 0 |
| AR | 117794 | 1 | 1123 | 31 | 176 |
| AS | 48271 | 0 | 2725 | 0 | 0 |
| CB | 28311 | 1 | 651 | 0 | 27 |
| CE | 2196 | 2 | 3553 | 8 | 0 |
| CL | 162015 | 692 | 60910 | 242 | 1 |
| CM | 148401 | 2593 | 35190 | 208 | 16 |
| CN | 42586 | 0 | 10310 | 0 | 49 |
| CT | 545352 | 1 | 4 | 3526 | 38731 |
| EX | 50164 | 249 | 23591 | 15 | 56 |
| GA | 111953 | 64 | 10391 | 82 | 0 |
| IB | 50393 | 154 | 8848 | 0 | 1 |
| MC | 78451 | 193 | 32224 | 66 | 90 |
| MD | 476969 | 22 | 210339 | 1 | 0 |
| ML | 4934 | 0 | 3816 | 0 | 0 |
| NC | 49894 | 64 | 10260 | 0 | 3 |
| PV | 133486 | 0 | 54594 | 0 | 142 |
| RI | 27307 | 59 | 2431 | 12 | 17 |
| VC | 314633 | 355 | 75799 | 0 | 182 |
De la tabla anterior llama mucho la atención el caso de Cataluña debido al gran número de desconocidos y al casi nulo número de test de antígenos que constrasta con el resto de CCAA. Es posible que pueda deberse a algún tipo de error en el volcado de los datos y que el número de desconocidos realmente se trate del número de diagnosticados con test de antígenos. Sin embargo, como no podemos constatarlo, no realizaremos cambio alguno. Lo que si haremos es eliminar la variable desconocido, ya que realmente son valores pérdidos con una importancia muy reducida sobre el total (excepto en el caso de Cataluña).
tipo_prueba<-select(tipo_prueba,-desconocido)Realizamos un análisis de correspondencias simple:
modelo_corresp<-CA(tipo_prueba,graph = FALSE)
print(modelo_corresp)## **Results of the Correspondence Analysis (CA)**
## The row variable has 19 categories; the column variable has 4 categories
## The chi square of independence between the two variables is equal to 355497.8 (p-value = 0 ).
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$col" "results for the columns"
## 3 "$col$coord" "coord. for the columns"
## 4 "$col$cos2" "cos2 for the columns"
## 5 "$col$contrib" "contributions of the columns"
## 6 "$row" "results for the rows"
## 7 "$row$coord" "coord. for the rows"
## 8 "$row$cos2" "cos2 for the rows"
## 9 "$row$contrib" "contributions of the rows"
## 10 "$call" "summary called parameters"
## 11 "$call$marge.col" "weights of the columns"
## 12 "$call$marge.row" "weights of the rows"Se rechaza la hipótesis nula de de independencia entre las CCAA y las diferentes pruebas diagnósticas (p-value=0), lo que indica que algunas CCAA difieren en las proporciones de casos diagnósticados por cada tipo de prueba respecto del total de casos.
# Valores propios
formattable(as.data.frame(round(get_eigenvalue(modelo_corresp),2)))| eigenvalue | variance.percent | cumulative.variance.percent | |
|---|---|---|---|
| Dim.1 | 0.09 | 91.10 | 91.10 |
| Dim.2 | 0.01 | 7.02 | 98.12 |
| Dim.3 | 0.00 | 1.88 | 100.00 |
fviz_screeplot(modelo_corresp,addlabels=TRUE)
Con las primeras dos dimensiones ya conseguimos explicar casi la totalidad (98,1%) de la variabilidad total de los datos. originales.
Vamos a representar gráficamente las constribuciones totales y relativas (cos2), tanto de los puntos fila (CCAA) como de los puntos columna (tipos de prueba), lo que nos indicará la importancia de cada punto en la construcción de cada eje y de la calidad de su representación de dichos puntos en el mapa.
# Análisis de contribuciones totales (importancia en la construcción de los ejes)
g31<-fviz_contrib(modelo_corresp,choice = "col",axes=1)+labs(title="Cont. test a Dim1",y="Contrib")
g32<-fviz_contrib(modelo_corresp,choice = "col",axes=2)+labs(title="Cont. test a Dim2",y="Contrib")
g33<-fviz_contrib(modelo_corresp,choice = "col",axes = 3)+labs(title="Cont. test a Dim3",y="Contrib")
g34<-fviz_contrib(modelo_corresp,choice = "row",axes = 1)+labs(title="Cont. CCAA a Dim1",y="Contrib")
g35<-fviz_contrib(modelo_corresp,choice = "row",axes = 2)+labs(title="Cont. CCAA a Dim2",y="Contrib")
g36<-fviz_contrib(modelo_corresp,choice = "row",axes = 3)+labs(title="Cont. CCAA a Dim3",y="Contrib")
(g31/g32/g33)|(g34/g35/g36)
# Análisis de contribuciones relativas (calidad de la representación)
g37<-fviz_cos2(modelo_corresp,choice = "col",axes=1)+labs(title="Cos2. test a Dim1",y="Cos2")
g38<-fviz_cos2(modelo_corresp,choice = "col",axes = 2)+labs(title="Cos2. test a Dim2",y="Cos2")
g39<-fviz_cos2(modelo_corresp,choice = "col",axes = 3)+labs(title="Cos2. test a Dim3",y="Cos2")
g310<-fviz_cos2(modelo_corresp,choice = "row",axes = 1)+labs(title="Cos2. CCAA a Dim1",y="Cos2")
g311<-fviz_cos2(modelo_corresp,choice = "row",axes = 2)+labs(title="Cos2. CCAA a Dim2",y="Cos2")
g312<-fviz_cos2(modelo_corresp,choice = "row",axes = 3)+labs(title="Cos2. CCAA a Dim3",y="Cos2")
(g37/g38/g39)|(g310/g311/g312)
Respecto a las contribuciones totales, Cataluña (respecto a los puntos fila) y el test de antígenos (respecto a los puntos columna) son los que más contribuyen en la construcción del primer eje. Esto se debe a la enorme diferencia que tiene Cataluña con las otras CCAA respecto a este tipo de test, que como hemos argumentado antes seguramente se deba a algún tipo de error. También contribuyen de forma significativa Madrid, Andalucía y Aragón, debido al escaso número de disgnosticados con test de antígenos (en términos relativos) frente al resto de CCAA. En cuanto a la segunda dimensión, se construye fundamentalmente a partir de Castilla la Mancha y el test de anticuerpos, debido al gran número de diagnosticados por esta prueba en esta comunidad frente al resto de CCAA. Por tanto, esta dimensión básicamente diferenciará a Castilla la Mancha y la prueba de anticuerpos del resto. Por último, la tercera dimensión se fundamentalmente por Cataluña y Aragón y por la prueba elisa. También contribuyen de forma importante Valencia, Galicia, Asturias y Cantabria.
Respecto a las contribuciones relativas, con las dos primeras dimensiones quedarán bien representados casi la totalidad de los puntos, a excepción de la prueba elisa y las comunidades de Valencia y Canarias (y en menor medida Baleares).
Realizamos 3 biplots: uno para CCAA para buscar similitudes entre las mismas, otro para las pruebas diagnósticas con el mismo objetivo y otro biplot simétrico para buscar relaciones entre CCAA y pruebas.
g313<-fviz_ca_row(modelo_corresp,repel=TRUE,col.row = "cos2",gradient.cols=c("#00AFBB","#E7B800","#FC4E07"))
g314<-fviz_ca_col(modelo_corresp,repel=TRUE,col.col = "cos2",gradient.cols=c("#00AFBB","#E7B800","#FC4E07"))
g315<-fviz_ca_biplot(modelo_corresp,arrow=c(TRUE,TRUE),repel=TRUE)
(g313|g314)/g315
A la vista de los gráficos anteriores, podemos concluir que existen CCAA que han apostado más que las otras por la prueba de antígenos, como Madrid,Castilla-León, País Vasco, Extremadura, Murcia o Andalucía, y muy especialmente Ceuta y Melilla. Por otro lado, Castilla la Mancha se diferencia del resto por el gran número de pruebas de anticuerpos que ha realizado. El resto se caracteriza por una preferencia mayor por la prueba pcr. Por último, la representación de Valencia, Canarias y la prueba elisa ya hemos explicado que no es fiable, ya que estan explicadas fundamentalmente por la tercera dimensión.
Hay que apuntar que, en términos absolutos, todas las CCAA (excepto Ceuta) han apostado de forma aplastante por la prueba pcr y que cuando se habla de preferencia de una comunidad, por ejemplo, por la prueba de antígenos quiere decir que la proporción sobre el total de casos de dicha comunidad es mayor que el de otra, y no que sea el tipo de prueba que más se ha hecho.
EJERCICIO 4
behavior<- read.spss("behavior.sav",to.data.frame = TRUE)
str(behavior)## 'data.frame': 15 obs. of 16 variables:
## $ IDFILA : Factor w/ 15 levels "Clase","Cita",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Correr : num 6.48 4 7.56 6.44 1.06 7.62 7.06 3.42 6.54 7.04 ...
## $ Hablar : num 2.79 0.44 0.92 0.48 0.58 5.71 0.54 0.81 4.02 0.75 ...
## $ Besar : num 6.9 0.27 4.73 4.08 1.29 6.62 7.92 4.25 2.79 3.83 ...
## $ Escribir : num 0.83 5.38 4.13 6.42 2 6.15 4.15 5.62 6.27 3.62 ...
## $ Comer : num 4.77 1.21 3.52 0.56 0.87 7.62 7.27 4.17 1.52 1.33 ...
## $ Dormir : num 5.4 5.23 1.96 6.71 3.37 7.23 8.25 7.54 4.92 6.1 ...
## $ Mascullar: num 5.38 5.88 3.83 6.46 3.6 5.48 7.69 4.04 4.87 2.79 ...
## $ Leer : num 1.73 6.12 1.83 5.04 1.23 5.42 6.52 4.19 7.27 4.29 ...
## $ Pelear : num 7.79 5.42 7.48 7.33 5.94 8.38 7.96 7.54 7.63 7.1 ...
## $ Eructar : num 7.23 6.77 6.85 6.5 4 7.58 7.79 6.19 6.42 3.96 ...
## $ Discutir : num 3.67 4.5 4.83 5.75 3.94 7.08 7.17 4.92 7.29 4.69 ...
## $ Saltar : num 7.21 4.58 5.88 6.71 1.58 7.29 7.52 5.46 6.69 5.25 ...
## $ Llorar : num 6.79 5.96 5.92 5.79 3.79 5.87 7.63 5.29 1.85 5.56 ...
## $ Reir : num 2.77 1 1.9 1.87 0.9 6.4 3.12 1.6 1.06 0.77 ...
## $ Gritar : num 7.06 5.21 6 7.04 2.08 7.67 7.35 4.12 6.58 4.87 ...
## - attr(*, "variable.labels")= Named chr [1:16] "" "" "" "" ...
## ..- attr(*, "names")= chr [1:16] "IDFILA" "Correr" "Hablar" "Besar" ...
## - attr(*, "codepage")= int 1252La muestra no cumple con las condiciones necesarias para realizar un AF, ya que se exige como mínimo una muestra de 50 observaciones. Incluso hay quien afirma que se necesitan aproximadamente 300 para obtener resultados útiles y relativamente estables (Tabachnick y Fidell, 2001). Por otro lado, se debería contar idealmente con 10 participantes por variable, y como mínimo con cinco por ítem (Nunnally y Bernstein, 1995).
Además, el hecho de contar con el mismo número de variables que de observaciones va a provocar que la matriz no sea invertible, al tener un determinante 0, lo que tampoco permitiría realizar el análisis.
behavior<-column_to_rownames(behavior,var="IDFILA")
cor_matrix<-cor(behavior) # Matriz de correlaciones de las columnas (comportamientos)
paste("Determinante: ",det(cor_matrix))## [1] "Determinante: 7.32583612683368e-26"cortest.bartlett(cor_matrix,n = 15)## $chisq
## [1] 472.6524
##
## $p.value
## [1] 4.562423e-48
##
## $df
## [1] 105KMO(cor_matrix)## Error in solve.default(r) :
## sistema es computacionalmente singular: número de condición recíproco = 8.96197e-18## matrix is not invertible, image not found## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor_matrix)
## Overall MSA = 0.5
## MSA for each item =
## Correr Hablar Besar Escribir Comer Dormir Mascullar Leer
## 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
## Pelear Eructar Discutir Saltar Llorar Reir Gritar
## 0.5 0.5 0.5 0.5 0.5 0.5 0.5Como se puede observar, el determinante de la matriz de correlaciones es cero, lo que tampoco permite calcular el índice KMO para comprobar la idoneidad del análisis factorial. Pese a todos los inconvenientes, vamos a eliminar una variable que esté muy correlacionada con otra y vamos a realizar el análisis.
corrplot(cor_matrix,method = "number",number.cex = 0.6)
Eliminamos la variable Gritar, que tiene una correlación de 0,96 con la variable Saltar. Calculamos determinante, test de Bartlett e índice KMO.
behavior2<-select(behavior,-Gritar)
cor_matrix2<-cor(behavior2)
paste("Determinante: ",det(cor_matrix2))## [1] "Determinante: 1.64721872449581e-10"cortest.bartlett(cor_matrix2,n = 15)## $chisq
## [1] 191.4775
##
## $p.value
## [1] 3.928727e-09
##
## $df
## [1] 91KMO(cor_matrix2)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor_matrix2)
## Overall MSA = 0.4
## MSA for each item =
## Correr Hablar Besar Escribir Comer Dormir Mascullar Leer
## 0.34 0.30 0.37 0.46 0.33 0.41 0.34 0.52
## Pelear Eructar Discutir Saltar Llorar Reir
## 0.42 0.54 0.41 0.43 0.31 0.77El determinante sigue siendo prácticamente cero, pero esta vez si es posible calcular el índice KMO. Aún así. nos da un valor muy pobre de 0,4 que indica que el análisis factorial no es adecuado en este caso. En contra, el test de Bartlett nos indica que sí podría realizarse el análisis.
Estimamos un modelo de tres factores por mínimos cuadrados generalizados ponderados con rotación varimax:
factorial<-fa(cor_matrix2,nfactors = 3,n.obs = 15,fm="gls",rotate="varimax")
factorial$loadings##
## Loadings:
## GLS1 GLS2 GLS3
## Correr 0.603 0.344 0.340
## Hablar 0.720 0.275 -0.161
## Besar 0.738 0.556
## Escribir 0.162 0.875
## Comer 0.779 0.349
## Dormir 0.265 0.608 0.526
## Mascullar 0.176 0.468 0.694
## Leer 0.897 0.201
## Pelear 0.717 0.365 0.435
## Eructar 0.343 0.423 0.677
## Discutir 0.472 0.785 0.268
## Saltar 0.605 0.390 0.499
## Llorar 0.192 0.864
## Reir 0.854 0.127 0.221
##
## GLS1 GLS2 GLS3
## SS loadings 4.154 3.459 3.142
## Proportion Var 0.297 0.247 0.224
## Cumulative Var 0.297 0.544 0.768Con tres dimensiones conseguimos explicar un 76,8% de las variables originales, lo cual es un valor aceptable.
A la vista de las saturaciones de la matriz rotada podemos agrupar los comportamientos en las dimensiones de la siguiente forma:
- Dimensión 1: Correr, Hablar, Besar, Comer, Pelear, Saltar y Reir.
- Dimensión 2: Escribir, Leer, Dormir y Discutir.
- Dimensión 3: Mascullar, Eructar y Llorar.
En términos generales y con alguna excepción podemos interpretar que la dimendión 1 recoge comportamientos más activos que en ciertos lugares podrían molestar a otros. La dimensión dos recoge comportamientos más tranquilos (a excepción de discutir) que no molestan a otros. Y la dimensión 3 recogería comportamientos más desagradables que se hacen de forma más íntima, como eructar o llorar.
Una vez interpretadas las dimensiones, vamos a ver como se comportan en cada uno de los lugares, mediante sus puntuaciones estimadas por el método de Thurstone:
factorial.puntuaciones<-factor.scores(behavior2,factorial,method = "Thurstone")
formattable(as.data.frame(factorial.puntuaciones$scores))| GLS1 | GLS2 | GLS3 | |
|---|---|---|---|
| Clase | 0.45358001 | -1.6556555 | 1.84128113 |
| Cita | -2.03491263 | 0.7651767 | 0.55775980 |
| Autobús | 0.47427080 | -0.6127469 | 0.44437825 |
| Cena | -0.36858610 | 0.7484725 | 0.63924803 |
| Parque | -1.05099979 | -0.5097940 | -0.06700956 |
| Iglesia | 2.64196012 | 0.6538362 | -0.54933377 |
| Entrevista | 0.74284076 | 0.5178370 | 1.64589495 |
| Acera | 0.11687315 | 0.2455945 | 0.02870380 |
| Cine | 0.13715288 | 1.8486048 | -1.26647882 |
| Bar | -0.04182687 | 0.2429192 | -0.22849666 |
| Ascensor | -0.10419808 | 0.7457495 | 0.59034890 |
| Baño | 0.38842215 | -0.1469316 | -0.47216534 |
| Dormitorio | -0.16884820 | -1.5673134 | -1.94538325 |
| Salón | -0.28003341 | -1.3228158 | -0.95809374 |
| Juego | -0.90569478 | 0.0470670 | -0.26065370 |
Los comportamientos englobados en la dimensión 1 serían más adecuados en una cita, en el parque o en un juego y no serían nada adecuados en la iglesia. Tiene sentido porque, por ejemplo, en una cita se come, habla, ríe o besa (o incluso pelea) y en el parque se corre o salta, comportamientos que no casan tan bien en una iglesia.
Por su parte, los comportamientos de la dimensión 2 (leer, escribir, dormir y discutir) son más adecuados en el clase, en el dormitorio o el salón y no son nada adecuados en el cine y en menor medida en la cena y el ascensor.
Por último, los de la tercera dimensión (eructar,llorar y mascullar) son más adecuados en el dormitorio, en el cine y el salón y no son nada adecuados en clase o en una entrevista.
EJERCICIO 5
breakfast<- read.spss("breakfast_overall.sav",to.data.frame = TRUE)
str(breakfast)## 'data.frame': 42 obs. of 16 variables:
## $ género: Factor w/ 2 levels "Hombre","Mujer": 1 2 1 2 1 2 1 2 1 2 ...
## $ TS : num 13 15 15 6 15 9 9 15 15 15 ...
## $ TM : num 12 11 10 14 9 11 14 10 12 13 ...
## $ MM : num 7 6 12 11 6 14 5 12 2 10 ...
## $ DJ : num 3 3 14 3 14 4 6 6 4 7 ...
## $ TC : num 5 10 3 7 13 7 8 9 5 6 ...
## $ BAM : num 4 5 2 8 2 6 4 2 8 4 ...
## $ RM : num 8 14 9 12 12 15 13 13 10 9 ...
## $ TMd : num 11 8 8 10 8 10 11 8 11 12 ...
## $ TMJ : num 10 9 7 9 7 8 12 7 3 11 ...
## $ TMg : num 15 12 11 15 10 12 15 11 13 14 ...
## $ PC : num 2 7 1 4 11 5 7 3 7 5 ...
## $ PD : num 1 1 6 1 1 2 2 1 9 2 ...
## $ DG : num 6 4 4 2 4 3 1 5 6 8 ...
## $ TCf : num 9 2 5 5 3 1 3 4 1 1 ...
## $ BM : num 14 13 13 13 5 13 10 14 14 3 ...
## - attr(*, "variable.labels")= Named chr [1:16] "Género" "Tostada sóla" "Tostada con mantequilla" "Magdalena y margarina" ...
## ..- attr(*, "names")= chr [1:16] "género" "TS" "TM" "MM" ...
## - attr(*, "codepage")= int 1252breakfast<-rename(breakfast,genero=género) # Quitamos la tilde de "género" para evitar problemasPara discriminar entre géneros según las preferencias de los individuos vamos a realizar un análisis discriminante lineal. Antes de proceder a estimar el modelo, realizaremos un pequeño análisis descriptivo mediante gráficos de densidad entre cada par de variables (alimentos) distinguiendo ambos géneros:
g51<-ggplot(breakfast)+geom_density(aes(TS,fill=genero),alpha=0.6)
g52<-ggplot(breakfast)+geom_density(aes(TM,fill=genero),alpha=0.6)
g53<-ggplot(breakfast)+geom_density(aes(MM,fill=genero),alpha=0.6)
g54<-ggplot(breakfast)+geom_density(aes(DJ,fill=genero),alpha=0.6)
g55<-ggplot(breakfast)+geom_density(aes(TC,fill=genero),alpha=0.6)
g56<-ggplot(breakfast)+geom_density(aes(BAM,fill=genero),alpha=0.6)
g57<-ggplot(breakfast)+geom_density(aes(RM,fill=genero),alpha=0.6)
g58<-ggplot(breakfast)+geom_density(aes(TMd,fill=genero),alpha=0.6)
g59<-ggplot(breakfast)+geom_density(aes(TMJ,fill=genero),alpha=0.6)
g510<-ggplot(breakfast)+geom_density(aes(TMg,fill=genero),alpha=0.6)
g511<-ggplot(breakfast)+geom_density(aes(PC,fill=genero),alpha=0.6)
g512<-ggplot(breakfast)+geom_density(aes(PD,fill=genero),alpha=0.6)
g513<-ggplot(breakfast)+geom_density(aes(DG,fill=genero),alpha=0.6)
g514<-ggplot(breakfast)+geom_density(aes(TCf,fill=genero),alpha=0.6)
g515<-ggplot(breakfast)+geom_density(aes(BM,fill=genero),alpha=0.6)
(g51|g52|g53)/(g54|g55|g56)/(g57|g58|g59)/(g510|g511|g512)/(g513|g514|g515)+ plot_layout(guides = 'collect')
En términos generales, existe un alto grado de solapamiento entre ambas distribuciones que ya nos indican que el análisis disciminante lineal podría no ser adecuado. Por ejemplo, la preferencia de ambos géneros respecto el alimento DJ es prácticamente idéntico, por lo que quizás sería adecuado eliminarlo del modelo. Situaciones parecidas se dan con TS y DG.
Por otro lado, hay variables que si podrían tener cierto poder discriminante. Por ejemplo, TM y TMg, donde hay una mayor preferencia por hombres que por mujeres; o PC y TCf, hacia donde hay una mayor preferencia por mujeres que por hombres.
Estimamos un primer modelo con todas las variables:
andisc<-lda(genero~.,data=breakfast)## Warning in lda.default(x, grouping, ...): variables are collinearandisc## Call:
## lda(genero ~ ., data = breakfast)
##
## Prior probabilities of groups:
## Hombre Mujer
## 0.5 0.5
##
## Group means:
## TS TM MM DJ TC BAM RM TMd
## Hombre 11.33333 6.571429 6.666667 8.47619 8.523810 6.666667 8.761905 8.47619
## Mujer 12.09524 9.285714 8.095238 8.47619 7.952381 5.904762 10.000000 8.52381
## TMJ TMg PC PD DG TCf BM
## Hombre 6.857143 9.380952 7.619048 5.285714 7.952381 6.571429 10.857143
## Mujer 8.095238 10.619048 6.238095 4.571429 7.142857 3.809524 9.190476
##
## Coefficients of linear discriminants:
## LD1
## TS -0.06241190
## TM 0.11206136
## MM 0.17207204
## DJ 0.02463572
## TC -0.02848316
## BAM 0.07466912
## RM -0.04674359
## TMd -0.04108829
## TMJ 0.08585917
## TMg -0.07464274
## PC -0.00483544
## PD 0.13975060
## DG -0.01269307
## TCf -0.28635549
## BM -0.12223426Al estimar el modelo, aparece un warning indicando que existe una alta colinealidad entre variables.
Para evaluar de forma correcta el poder predictivo del modelo sería necesario contar con un conjunto test de datos, sin embargo, sólo contamos con 42 observaciones, por lo que la evaluación se realizará sobre los mismos datos de entrenamiento.
probs<-predict(andisc,breakfast,type="prob") # Estimación de probabilidades con el modelo realizado
table(probs$class,breakfast$genero) # Tabla de confusión##
## Hombre Mujer
## Hombre 17 7
## Mujer 4 14mean(probs$class==breakfast$genero) # Accuracy## [1] 0.7380952Obtenemos una precisión algo pobre del 73,8%. Sin embargo, mejora la precisión a priori, es decir, considerendo todas las observaciones pertenecientes a la clase mayoritaria (50%).
Analizamos las correlaciones:
correlaciones<-cor(breakfast[,-1])
corrplot(correlaciones,method = "number",type ="upper",diag = FALSE,number.cex =0.7)
Existen variables con correlaciones algo altas, como TM, DJ o RM.
Vamos a estimar un segundo modelo en el que eliminamos las variables TS,DJ y DG, que como hemos visto en los gráficos de densidad apenas tienen poder discriminante y lo único que pueden estar haciendo es introducir ruido en el modelo.
breakfast2<-breakfast %>% dplyr::select(-DJ,-TS,-DG) # Eliminamos variables TS,DJ,DG
andisc2<-lda(genero~.,data=breakfast2)
probs2<-predict(andisc2,breakfast2,type="prob")
table(probs2$class,breakfast2$genero)##
## Hombre Mujer
## Hombre 17 7
## Mujer 4 14mean(probs2$class==breakfast2$genero)## [1] 0.7380952Obtenemos la misma precisión que el modelo anterior, por lo que en base al principio de parsimonia nos quedamos con el segundo modelo que es más simple.
Además, con este modelo no nos aparece el warning por alta colinealidad.
En definitiva, las preferencias en el desayuno parecen tener cierto poder discriminante entre hombres y mujeres.