Enunciado

El objetivo de esta práctica es obtener las componentes principales para la variable ozono a partir de la combinación de la temperatura, la radiación social y la velocidad del viento. Para ello, dispones de una tabla de datos que recoge las lecturas de estas variables para la estación meteorológica de Valencia-viveros (Estación: 46250043-València - Vivers).

La tabla de datos es una serie temporal del 01/01/2010 a 31/12/2010. Los datos son de frecuencia diaria.

Los datos disponibles en la tabla son los siguientes:

Tipo de Variable Descripcion
Ozono Contínua Nivel medio diario de O3 en µg/m³
RSolar Contínua Radiación ultravioleta efectiva. Promedio diario entre las 8:00 y las 11:59. Medida en mW/m2
Temp Contínua Temperatura media mensual. Grados centígrados
Veloc Contínua Velocidad media diario del viento Km/h

Se pide: 1. Calcula la matriz de correlaciones, represéntala gráficamente. 2. A partir de la anterior calcula los valores y vectores propios. 3. Obtén las componentes resultantes y represéntalas gráficamente. 4. Finalmente, representa la incidencia de cada variable en las componentes obtenidas. 5. Trata de interpretar las componentes que has obtenido.

Resolución de la práctica


1. Calcula la matriz de correlaciones, represéntala gráficamente.

En primer lugar se cargan los datos de la ubicación disponible en la web del grupo IMEQ

url <- "http://www.uv.es/imeq/files/ozovlc.txt"
ozovlc <- read.table(url, header = T)
ozovlc <- ozovlc [,c(3,1,2,4)]
names (ozovlc) <- c ("Ozono","Temp","Veloc","RSolar")
library (corrplot)
S<-cor(ozovlc, use="complete.obs")
corrplot.mixed(S)

2. A partir de la anterior calcula los valores y vectores propios.

rdo<-eigen(S)

#Valores propios
round(rdo$values,3)
## [1] 2.245 1.069 0.515 0.171
#Vectores propios
vp <- data.frame(rdo$vectors)         
row.names(vp) <- names (ozovlc)
names(vp) <- paste0("Comp.",c(1:4))
round(vp,3)
##        Comp.1 Comp.2 Comp.3 Comp.4
## Ozono  -0.522 -0.277  0.754 -0.287
## Temp   -0.570  0.261 -0.520 -0.580
## Veloc  -0.157 -0.898 -0.397  0.107
## RSolar -0.615  0.222 -0.056  0.755

3. Obtén las componentes resultantes y represéntalas gráficamente.

A partir de la matriz de correlación obtenida en el apartado anterior, obtenemos las componentes principales:

Resumen con % varianza explicada

CP<-princomp(covmat=S)
summary(CP)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3    Comp.4
## Standard deviation     1.4984973 1.0338308 0.7174819 0.4134242
## Proportion of Variance 0.5613735 0.2672015 0.1286951 0.0427299
## Cumulative Proportion  0.5613735 0.8285750 0.9572701 1.0000000

Si queremos ver la estructura de la salida

str(CP)
## List of 7
##  $ sdev    : Named num [1:4] 1.498 1.034 0.717 0.413
##   ..- attr(*, "names")= chr [1:4] "Comp.1" "Comp.2" "Comp.3" "Comp.4"
##  $ loadings: loadings [1:4, 1:4] -0.522 -0.57 -0.157 -0.615 -0.277 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "Ozono" "Temp" "Veloc" "RSolar"
##   .. ..$ : chr [1:4] "Comp.1" "Comp.2" "Comp.3" "Comp.4"
##  $ center  : num [1:4] NA NA NA NA
##  $ scale   : Named num [1:4] 1 1 1 1
##   ..- attr(*, "names")= chr [1:4] "Ozono" "Temp" "Veloc" "RSolar"
##  $ n.obs   : logi NA
##  $ scores  : NULL
##  $ call    : language princomp(covmat = S)
##  - attr(*, "class")= chr "princomp"

Varianzas de las CP (valores propios)

CP$sdev**2
##    Comp.1    Comp.2    Comp.3    Comp.4 
## 2.2454940 1.0688061 0.5147803 0.1709196

Pesos para formar las CP, (loadings)

print(loadings(CP),cutoff=0)
## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4
## Ozono  -0.522 -0.277  0.754 -0.287
## Temp   -0.570  0.261 -0.520 -0.580
## Veloc  -0.157 -0.898 -0.397  0.107
## RSolar -0.615  0.222 -0.056  0.755
## 
##                Comp.1 Comp.2 Comp.3 Comp.4
## SS loadings      1.00   1.00   1.00   1.00
## Proportion Var   0.25   0.25   0.25   0.25
## Cumulative Var   0.25   0.50   0.75   1.00

Observamos que con las tres primeras dos primeras componentes aportan más del 82% de la varianza, mientras que con la tercera se aumenta un 13% hasta el 95%. La última componente agrega un escaso 4.3%.

Gráficamente puede representarse el peso de las componentes principales por pares:

G<-CP$loadings
plot(G[,1:2],type="n"
     ,xlim=c(-0.7,0.3), cex.main=.7, cex.lab =.5
     )
for (i in 1:nrow(G)) arrows(0,0,G[i,1],G[i,2])
text(G[,1:2],dimnames(G)[[1]],cex=.9,adj=0)
title("Pesos de las C.P. 1 y 2")

plot(G[,2:3],type="n",xlim=c(-1,0.4) , cex.main=.7, cex.lab =.5)
for (i in 1:nrow(G)) arrows(0,0,G[i,2],G[i,3])
text(G[,2:3],dimnames(G)[[1]],cex=.9,adj=0)
title("Pesos de las C.P. 2 y 3")

Vemos que la componentes radiación solar y temperatura apuntan con intensidad parecida hacia la misma dirección por lo que podrían estar aportando la misma información.

4. Finalmente, representa la incidencia de cada variable en las componentes obtenidas.

Gráficamente observamos la incidencia de cada variable en cada una de las cuatro componentes:

barplot(CP$loadings[,1],main="Componente 1",beside=T,col="black")
barplot(CP$loadings[,2],main="Componente 2",beside=TRUE,col="grey")
barplot(CP$loadings[,3],main="Componente 3",beside=TRUE,col="lightgrey")
barplot(CP$loadings[,4],main="Componente 4",beside=TRUE,col="orange")

5. Trata de interpretar las componentes que has obtenido.

La primera componente (con una capacidad explicativa del 56%) identifica días donde radiación solar, ozono y temperatura muestran valores altos siendo, a su vez, días poco ventosos. La segunda componente (con una capacidad explicativa del 27%) estaría identificando días ventosos, dado que el resto de variables no toman valores significativamente altos. En la tercera componente, estaríamos identificando días en los que el nivel de ozono toma valores elevados, mientras que la temperatura y el viento influyen negativamente. La cuarta componente, por su bajo aporte a la variabilidad del modelo, la podríamos dejar fuera, reduciendo el número de variables a considerar a tres.

La decisión de elegir las tres primeras componentes, se corrobora con son días, a su vez, por el gráfico del bastón roto, donde observamos un cambio de pendiente en la tercera variable, indicando que la variabilidad explicada por la cuarta variable es baja (un 4,5%).

plot(CP,type="l",main="Valores propios")