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.
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)
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
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.
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")
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")