Análisis de las emisiones de CO2 a nivel mundial
Jose Luis Canillas Cambronero
20 de noviembre de 2018
1. Introducción.
El objetivo del trabajo va a ser estudiar a nivel mundial y europeo, la evolución de una serie de variables relacionadas con el medio ambiente y la energía, centrando sobretodo el análisis en las emisiones de gases de CO2 responsables del cambio climático.
Concretando más, el trabajo va a consistir en tres partes diferenciadas, aunque relacionadas:
En primer lugar un análisis de la varianza referido a la relación entre el PIB per cápita y el precio de los combustibles fósiles.
Un análisis de componentes principales, en el cual vamos a intentar mostrar la interrelación que existe entre variables de crecimiento, consumo, industria, con variables contaminantes.
Un análisis cluster, en el cual se va a estudiar, para el caso de la UE, la evolución de la producción de renovables y emisiones de CO2 desde 1990 hasta 2016, realizando agrupaciones según como avancen estas dos variables.
Se pretende dar una visión general de lo importante que es tener en cuenta las emisiones de gases nocivos, y como la utilización de energías renovables es crucial para luchar contra el cambio climático (esto lo veremos muy bien en el análisis clúster).
2. Análisis de la varianza.
En primer lugar abrimos los datos y transformamos las variables nominales en factor.
library(readxl)
datos<- read_excel("D:/Universidad/4_ADE/Tecnicas_estadisticas/Trabajo_individual/datos_VAR.xlsx")
datos$P_GSO_N <- as.factor(datos$P_GSO_N)
datos$P_DIE_N <- as.factor(datos$P_DIE_N)2.1. Análisis exploratorio.
Resulta interesante realizar un análisis exploratorio de los datos para ver si alguna variable contiene algún error o algún valor perdido:
library(dplyr)
datos_su <- datos %>%
select(-c(Pais,P_GSO_N,P_DIE_N))
knitr::kable(summary(datos_su))| PIB_per_capita | P_gasolina | P_diesel | PIB_PC_LN | |
|---|---|---|---|---|
| Min. : 1813 | Min. :0.600 | Min. :0.250 | Min. : 7.503 | |
| 1st Qu.: 9227 | 1st Qu.:1.270 | 1st Qu.:1.090 | 1st Qu.: 9.130 | |
| Median :14119 | Median :1.570 | Median :1.510 | Median : 9.555 | |
| Mean :23694 | Mean :1.509 | Mean :1.384 | Mean : 9.689 | |
| 3rd Qu.:41375 | 3rd Qu.:1.800 | 3rd Qu.:1.660 | 3rd Qu.:10.630 | |
| Max. :89275 | Max. :2.270 | Max. :2.110 | Max. :11.399 |
Vemos como no existen valores perdidos en la base de datos. Además, el PIB per cápita se ha expresado en logaritmos para eliminar la influencia de posibles valores atípicos. También tenemos la variable precio de la gasolina y precio del diesel categorizadas en precio Alto y Bajo, el objetivo es desglosar aquellos paises cuyo coste de adquisición del combustible es mayor o menor.
El precio de los combustibles vienen expresados en dólares, al igual que el PIB per cápita, aunque en este último, al estar expresados en logaritmos, se van a tomar los valores más altos como referencia a los países con mayor nivel de PIB per cápita.
Se va a estudiar a continuación una posible relación entre el precio de la gasolina y el diésel (alto o bajo, considerando alto si supera los 1,5 dólares) y el PIB per cápita. Gráficamente lo podemos observar en las siguientes cajas:
boxplot(datos$PIB_PC_LN ~ datos$P_GSO_N,
main = "Precio de la gasolina según el PIB per cápita", las = 1, col = "deepskyblue1",
ylab = "PIB per cápita (Log)", xlab = "Precio gasolina (Nominal)")
boxplot(datos$PIB_PC_LN ~ datos$P_DIE_N,
main = "Precio del diésel según el PIB per cápita", las = 1, col = "deepskyblue1",
ylab = "PIB per cápita (Log)", xlab = "Precio gasolina (Nominal)")
A simple vista vemos como aquellos países que tengan un mayor nivel de renta per cápita, a su vez van a poseer mayores gastos en adquisición de combustibles. Existen algunas excepciones vistas en las colas del gráfico. También vemos como la relación en el caso de la gasolina y el diésel es muy similar, aún así vamos a incluir estas dos variables como independiente a la hora de estimar el modelo.
2.2. Estimación y validación del modelo.
El modelo estimado va a contrastar si realmente existen diferencias entre los grupos. La variable dependiente serán las emisiones de CO2 y las variables independientes el precio de la gasolina y el precio del diésel. Veremos si existen diferencias en la cantidad emitida según las energías renovables que se produzcan.
modelo1 <- aov(datos$PIB_PC_LN ~ datos$P_GSO_N, datos$P_DIE_N)
summary(modelo1)## Df Sum Sq Mean Sq F value Pr(>F)
## datos$P_GSO_N 1 12.30 12.301 19.45 5.99e-05 ***
## Residuals 47 29.73 0.633
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p-valor es inferior a 0,05 lo que en principio podemos rechazar la hipótesis nula de medias iguales y afirmar que si existen diferencias entre los diferentes grupos. El precio de la gasolina llega a incidir, al menos en parte, el nivel de renta per cápita. Realmente ya para validar lo anterior dicho se debe contrastar el test de la F mediante el cumplimiento de las hipótesis de homogeneidad de varianzas y normalidad de los residuos.
Previamente, podemos verlo o intuirlo gráficamente:
par(mfrow = c(2, 2))
plot(modelo1)
En el primer gráfico se observa cómo no existe mucha dispersión de los residuos, ya que no realiza un posible forma de embudo, lo que en cierto modo confirma la homogeneidad de las varianzas. En el segundo gráfico los puntos se acercan bastante a la línea, esto puede confirmar la normalidad de los residuos.
2.3. Validación de hipótesis.
Lo dicho anteriormente se puede contrastar con varias pruebas, en concreto vamos a determinar si verdaderamente existe homogeneidad en las varianzas y normalidad en los residuos.
bartlett.test(datos$PIB_PC_LN ~ datos$P_GSO_N, datos$P_DIE_N)##
## Bartlett test of homogeneity of variances
##
## data: datos$PIB_PC_LN by datos$P_GSO_N
## Bartlett's K-squared = 0.82779, df = 1, p-value = 0.3629shapiro.test(residuals(modelo1))##
## Shapiro-Wilk normality test
##
## data: residuals(modelo1)
## W = 0.98188, p-value = 0.6462Vemos como el test de Bartlett es superior a 0,05 lo que no llegamos a rechazar la hipótesis nula de homogeneidad en las varianzas. Justo lo que sospechábamos en el gráfico de arriba se cumple.
Por otro lado, el test de Shapiro tampoco rechaza la hipótesis nula de normalidad en los residuos, ya que es muy superior a 0,05 y por tanto afirmamos que los residuos están normalizados.
Por último, no sería necesario realizar el test de Kruskal ya que se han cumplido las hipótesis anteriores y aunque solo existe un grupo de niveles (Gasolina a precio alto o bajo), vamos a ver con el test de Tukey las diferencias existentes y si es significativa.
2.4. Relación entre PIB per cápita y precio de los combustibles.
Como se ha dicho anteriormente, procedemos a realizar el test de Tukey:
TukeyHSD(modelo1, conf.level = 0.95)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = datos$PIB_PC_LN ~ datos$P_GSO_N, data = datos$P_DIE_N)
##
## $`datos$P_GSO_N`
## diff lwr upr p adj
## Bajo-Alto -1.003947 -1.461948 -0.545946 5.99e-05Confirma todo lo anterior, ya que la relación es significativa y por tanto existen diferencias entre el grupo Alto y Bajo. Esto quiere decir que aquellos países que tengan una renta per cápita mayor, normalmente también van a tener precios de los combustibles muy elevados, lo mismo, pero a la inversa ocurre con aquellos países que tenga una renta per cápita menor, que tendrán precio de la gasolina menores. Esto puede influir en las emisiones o en la producción de energías no renovables, ya que, al tener un precio mayor, se optan por otros tipos de combustibles. Realmente esto no tiene por qué ser así y la producción o no de energías renovables puede verse condicionada por otro tipo de factores. Mas tarde y con el análisis de los componentes principales veremos que existe relación entre el PIB y el consumo de un país y sus emisiones de CO2.
3. Análisis de Componentes Principales.
En este apartado se va a intentar reducir la dimensión de un conjunto de variables para intentar dar una explicación sencilla (como mucho utilizar dos dimensiones) y así explicar que variables de las seleccionadas influyen en gran medida en las emisiones de CO2.
En este caso vamos a utilizar una base de datos distinta, la cual contiene un amplio conjunto de variables referidas tanto a emisiones, como variables de crecimiento económico e industria (que son las que utilizaremos en esta técnica).
rm(list=ls())
library(readxl)
datos2 <- read_excel("datos_ACP.xlsx")3.1. Análisis exploratorio.
A continuación se muestra un análisis exploratorio de los datos:
datos2_su <- datos2 %>% select(-Pais)
pander::pandoc.table(summary(datos2_su))| ECO2_LIQUIDO | VA_INDUS | GTO_CONSUFINAL | PIB |
|---|---|---|---|
| Min. : 792.1 | Min. :11.29 | Min. :3.417e+09 | Min. :3.840e+09 |
| 1st Qu.: 8221.4 | 1st Qu.:20.23 | 1st Qu.:3.309e+10 | 1st Qu.:5.786e+10 |
| Median : 31323.5 | Median :25.19 | Median :1.443e+11 | Median :2.518e+11 |
| Mean : 113092.1 | Mean :25.88 | Mean :5.841e+11 | Mean :9.635e+11 |
| 3rd Qu.: 89628.8 | 3rd Qu.:31.32 | 3rd Qu.:3.150e+11 | 3rd Qu.:5.356e+11 |
| Max. :2114139.2 | Max. :52.78 | Max. :1.105e+13 | Max. :1.621e+13 |
| CONTAMI_AIRE | EMI_CO2 |
|---|---|
| Min. : 5.862 | Min. : 1984 |
| 1st Qu.:12.334 | 1st Qu.: 30678 |
| Median :17.909 | Median : 64602 |
| Mean :19.794 | Mean : 289143 |
| 3rd Qu.:21.121 | 3rd Qu.: 248315 |
| Max. :99.343 | Max. :5254279 |
Vemos como por ejemplo en la variable “EMI_CO2” las observaciones son muy heterogéneas, y el máximo y el mínimo difieren muchísimo. Lo importante del análisis será la variabilidad que exista y no la escala mayor de estas, por ello se tipificarán los datos posteriormente.
De la base de datos vamos a escoger un subconjunto para el análisis de componentes, este estará formado por todas las variables de la base a excepción de “EMI_CO2”, que posteriormente se integrará como variable dependiente en un modelo de regresión.
subdatos <- subset(datos2, select = c(ECO2_LIQUIDO,
VA_INDUS,GTO_CONSUFINAL,
PIB,CONTAMI_AIRE))La variables referidas a emisiones de CO2 están expresadas en toneladas métricas, el VAB de la industria en porcentaje del PIB, las variable PIB y Gasto en consumo final están expresadas en dólares y la variable contaminación del aire es el porcentaje de personas que exceden el valor indicativo de la Organización Mundial de la Salud (OMS).
3.2. Análisis de las correlaciones.
Es importante que las variables originales estén correlacionadas entre sí, por ello lo primero que se va a realizar es la matriz de correlaciones:
correlaciones <- cor(subdatos)
knitr::kable(correlaciones)| ECO2_LIQUIDO | VA_INDUS | GTO_CONSUFINAL | PIB | CONTAMI_AIRE | |
|---|---|---|---|---|---|
| ECO2_LIQUIDO | 1.0000000 | -0.0651241 | 0.9817992 | 0.9730240 | -0.1372631 |
| VA_INDUS | -0.0651241 | 1.0000000 | -0.1222778 | -0.1205091 | 0.4952104 |
| GTO_CONSUFINAL | 0.9817992 | -0.1222778 | 1.0000000 | 0.9973084 | -0.1679659 |
| PIB | 0.9730240 | -0.1205091 | 0.9973084 | 1.0000000 | -0.1800074 |
| CONTAMI_AIRE | -0.1372631 | 0.4952104 | -0.1679659 | -0.1800074 | 1.0000000 |
Podemos encontrar algunas relaciones a simple vista como GTO_CONSUFINAL con ECO2_LIQUIDO o PIB, aunque de manera más visual lo vemos gráficamente a continuación:
library(corrplot)## corrplot 0.84 loadedsig <- cor.mtest(subdatos, conf.level = 0.95)
corrplot(correlaciones, method = "number", type = "lower",p.mat = sig[[1]])
Se muestran las variables menos correlacionadas tachadas con una cruz, dejando visuales aquellas variables que sí tienen más relación con otras variables. Encontramos por ejemplo que el PIB tiene una correlación prácticamente perfecta con la variable gasto en consumo final. También se observa relación, aunque menor entre el VAB de Industria y la contaminación del aire y sobre todo entre las emisiones de CO2 líquido y el gasto en consumo final y el PIB.
Esta visualización también la encontramos con la función ezCor:
ez::ezCor(subdatos, r_size_lims = c(4, 8), label_size = 3)
En la parte superior se muestran las correlaciones de forma numérica, y en la parte inferior se muestran los gráficos de dispersión con sus respectivos intervalos de confianza, así como una recta de regresión donde se acumulan nuestros casos (cuanto más próximos estén los puntos a la recta, mayor correlación).
Antes de entrar en la creación de los componentes principales, se va a proceder a tipificar los datos para evitar problemas de escala.
DatosTIP <- scale(subdatos)
knitr::kable(head(DatosTIP))| ECO2_LIQUIDO | VA_INDUS | GTO_CONSUFINAL | PIB | CONTAMI_AIRE |
|---|---|---|---|---|
| -0.0630390 | -0.2089810 | -0.1789591 | -0.2164126 | -0.4211889 |
| -0.3742421 | -0.0793207 | -0.3602447 | -0.3965252 | 0.3612120 |
| -0.2761161 | -0.0673243 | -0.2316788 | -0.2306296 | -0.2279063 |
| 0.7533817 | -0.7057383 | 0.5901459 | 0.6081278 | -0.5803797 |
| -0.3402006 | -0.3024668 | -0.3444102 | -0.3792009 | 0.5318108 |
| 0.4509090 | 0.2160894 | 0.2654298 | 0.3419087 | -0.8620341 |
3.3. Estimación del modelo de componentes principales.
Ahora si procedemos a reducir la dimensionalidad de las variables una vez que hemos visto que existe correlación entre estas.
modelo1 <- prcomp(subdatos, scale = TRUE)
summary(modelo1)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.7431 1.1955 0.7092 0.1673 0.03965
## Proportion of Variance 0.6077 0.2859 0.1006 0.0056 0.00031
## Cumulative Proportion 0.6077 0.8935 0.9941 0.9997 1.00000En el summary se muestran las cinco componentes objeto del análisis, cada una de ellas contiene tres tipos de información muy interesantes. Nos interesa sobre todo la proporción de varianza que explica cada componente, ya que esto nos hace una idea de cuantos componentes vamos a escoger. Vemos en la proporción acumulativa que si recogemos dos componentes (PC1 y PC2) ya estaríamos explicando un 89% del modelo.
Para una mayor exactitud, vamos a seleccionar a aquellas componentes cuyo autovalor sea mayor que 1.
autovalores <- modelo1$sdev^2
autovalores## [1] 3.038244040 1.429275006 0.502920685 0.027987922 0.001572348selec <- sum(autovalores > 1)
selec## [1] 2Hemos procedido buscando en el modelo aquellos autovalores mayores que 1, después seleccionamos aquellos que superen la unidad, y tal como sospechábamos anteriormente, nos quedamos con dos componentes. Gráficamente lo podemos observar en el gráfico de sedimentación:
plot(autovalores, main = "Gráfico de Sedimentación",
xlab = "Nº de Autovalor", ylab="Valor", pch = 16, col = "red4", type = "b",
lwd = 2, las = 1)
abline(h = 1, lty = 2, col = "green4")
3.4. Cargas factoriales.
Una vez seleccionadas las componentes principales, vamos a medir la relación entre las componentes retenidas y las variables originales, a esto se le denomina “cargas factoriales”.
coeficientes <- modelo1$rotation
cargas <- t(coeficientes[, 1:2])*(sqrt(autovalores[1:2]))
knitr::kable(cargas)| ECO2_LIQUIDO | VA_INDUS | GTO_CONSUFINAL | PIB | CONTAMI_AIRE | |
|---|---|---|---|---|---|
| PC1 | 0.9750426 | -0.2188793 | 0.9900564 | 0.9882515 | -0.2877069 |
| PC2 | -0.1770500 | -0.8406113 | -0.1273915 | -0.1204845 | -0.8127459 |
Para el componente PC1 las mayores cargas se concentran en las variables ECO2_LIQUIDO, GTO_CONSUFINAL y PIB. Para el componente PC2 las mayores cargas se concentran en las variables VA_INDUS y CONTAMI_AIRE.
Gráficamente lo podemos observar mejor:
factoextra::fviz_pca_biplot(modelo1)
A partir de aquí, vamos a nombrar a cada componente en función de las cargas factoriales de cada una de estas, así pues, tenemos que:
- PC1 va a venir referido a “Factores económicos ligados a la contaminación”.
- PC2 va a venir referido a “Contaminación provocada por la industria”.
Por último, guardamos las dos componentes retenidas en la base de datos original para realizar a continuación la regresión lineal:
datos2 <- cbind(datos2, modelo1$x[, 1:2])3.5. Regresión lineal: comportamiento de las emisiones de CO2.
Por último, vamos a crear un regresión lineal donde la variable dependiente será “EMI_CO2”, que viene referida a las emisiones de CO2 y las variables independientes los componentes nombrados en el apartado anterior.
regre2 <- lm(EMI_CO2 ~ PC1+PC2, data = datos2)
summary(regre2)##
## Call:
## lm(formula = EMI_CO2 ~ PC1 + PC2, data = datos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -339890 -40425 -4653 34152 1039594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 289143 24650 11.730 5.73e-16 ***
## PC1 415281 14277 29.087 < 2e-16 ***
## PC2 -123827 20816 -5.949 2.63e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 179500 on 50 degrees of freedom
## Multiple R-squared: 0.9463, Adjusted R-squared: 0.9442
## F-statistic: 440.7 on 2 and 50 DF, p-value: < 2.2e-16Observamos como el resultado de la regresión es significativo, ya que hemos obtenidos un p-valor muy bajo y obtenemos por tanto un nivel de confianza del 99% para los dos componentes, además de un R2 muy alto. También observamos el valor de la estimación, el cual es positivo para los factores económicos ligados a la contaminación y negativo para la contaminación ligada a la industria.
Esto quiere decir en primer lugar, que las emisiones de CO2 van a estar ligadas de forma positiva con las variables dentro de los factores económicos, si se produce una expansión del PIB o del gasto en consumo final, también se producirá un aumento de las emisiones de CO2, especialmente las compuestas por aquellas derivadas de líquidos. En segundo lugar, tendríamos una relación negativa entre las emisiones de CO2 y contaminación ligada a la industria, en términos de que, si el valor agregado bruto de la industria es menor, también se produce una reducción de las emisiones de CO2 y también se produce dicha reducción en la cantidad de contaminación en el aire.
Vemos así una relación clara tanto en variables económicas como industriales, una posible causa de las emisiones de CO2 en nuestro planeta.
4. Análisis clúster.
Con el análisis clúster, en vez de tener países como observaciones, ahora tenemos años, y la evolución va a venir referida a la UE28. Vamos a agrupar la producción de renovables y emisiones de CO2 en diferentes grupos, según creamos conveniente y según lo que arrojen los análisis. Aclarar que las emisiones de CO2 vienen expresadas en toneladas de petróleo equivalente, al igual que la producción de renovables.
En primer lugar, como siempre, vamos a abrir los datos, además de colocar correctamente la primera columna, como los nombres de las observaciones:
library(readxl)
datos3 <- data.frame(read_excel("Datos_CLUS.xlsx"))
datos3$Fecha <- as.factor(datos3$Fecha)
rownames(datos3) <- datos3[,1]
datos3[,1] <- NULL4.1. Análisis exploratorio.
Se muestran a continuación las dos variables que vamos a utilizar para el análisis, asi como una serie de gráficos exploratorios que nos sirven para ponernos en situación.
pander::pandoc.table(summary(datos3))| Renovables | Emisiones_CO2 |
|---|---|
| Min. :38.95 | Min. :3614 |
| 1st Qu.:47.58 | 1st Qu.:4011 |
| Median :54.83 | Median :4304 |
| Mean :60.44 | Mean :4203 |
| 3rd Qu.:77.95 | 3rd Qu.:4406 |
| Max. :88.95 | Max. :4545 |
par(mfrow = c(2, 2))
hist(datos3$Renovables, main = "Producción de renovables", xlab = "")
hist(datos3$Emisiones_CO2, main = "Emisiones de CO2", xlab = "")
boxplot(datos3$Renovables)
boxplot(datos3$Emisiones_CO2)
En los gráficos observamos claramente como se sitúan las producciones de energía renovable y emisiones de CO2 a lo largo de los últimos 27 años. La mediana se sitúa a niveles muy bajos para el caso de producción de renovables, en su lugar, la mediana de las emisiones se mantiene muy alta en el tiempo. Si es verdad que el esfuerzo en contaminar menos en los últimos años es favorable, pero el resultado si tenemos en cuenta en periodo de 27 años es muy desfavorable.
La media de renovables se sitúa en 60,44 y la media de contaminación en 4203. Se ve un avance si tenemos en cuenta el mínimo de renovables y el máximo de emisiones.
4.2. Elección óptima de clústeres: Método jerárquico.
En primer lugar, vamos a utilizar el método de clasificación jerárquico. También usaremos la función NbClust para ver cuántos grupos sería óptimo formar. En combinación entre nuestro criterio y los distintos análisis que se desarrollen, se decidirá formar más o menos grupos.
Vamos a tipificar los datos para evitar problemas de escala:
datostip <- scale(datos3)
pander::pandoc.table(head(datostip))| Renovables | Emisiones_CO2 | |
|---|---|---|
| 1990 | -1.251 | 1.216 |
| 1991 | -1.158 | 1.008 |
| 1992 | -1.182 | 0.5133 |
| 1993 | -1.011 | 0.2279 |
| 1994 | -1.04 | 0.1667 |
| 1995 | -0.941 | 0.3495 |
Para crear el clúster para vamos a utilizar el método de Ward, para calcular las distancias vamos a utilizar la distancia euclidea.
Previamente podemos observar en el siguiente gráfico de dispersión como situan las observaciones:
library(scales)
par(mfrow = c(1,1 ))
plot(datostip, col = alpha("orange", 0.5), pch = 19, las=1,
main = "Relación Emisiones CO2 y Renovables por año")
text(datostip, rownames(datostip), pos = 3, cex = .6)
Existe una relación entre una mayor contaminación y una menor utilización de energias renovables. Asi como una relación entre una mayor utilización de energías renovables y una menor contaminación.
A continuación se muestra el dendrograma:
d <- dist(datostip, method = "euclidean")
cluster <- hclust(d, method="ward.D")
plot(cluster, cex = .6, xlab = "", ylab = "Distancia",
sub = "Cluster de grupos para los años") 
A priori podríamos podar el dendrograma en dos grupos, pero uno de ellos tendría excesivas observaciones comparado con el otro. Vamos a ver pues, que nos dice la función NBClust:
NbClust::NbClust(data = datostip,
distance = "euclidean",
method = "ward.D", max.nc = 5)
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 13 proposed 2 as the best number of clusters
## * 6 proposed 3 as the best number of clusters
## * 2 proposed 4 as the best number of clusters
## * 3 proposed 5 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 2
##
##
## *******************************************************************## $All.index
## KL CH Hartigan CCC Scott Marriot TrCovW TraceW Friedman
## 2 12.6385 105.1881 17.9502 2.3584 59.4584 95.6424 30.4748 9.9856 8.8409
## 3 1.0040 95.3588 13.4165 1.6114 98.6667 50.3694 19.3345 5.8123 24.3182
## 4 1.3608 99.2672 9.8343 2.4712 117.8209 44.0505 7.3513 3.7282 32.6874
## 5 0.7390 104.0146 10.6826 2.5756 133.8737 37.9805 2.8602 2.6115 42.5123
## Rubin Cindex DB Silhouette Duda Pseudot2 Beale Ratkowsky Ball
## 2 5.2075 0.3515 0.4436 0.7114 0.4121 24.2520 1.3473 0.6355 4.9928
## 3 8.9466 0.3402 0.6708 0.5624 0.2781 15.5760 2.2251 0.5440 1.9374
## 4 13.9479 0.4182 0.6522 0.5273 0.3877 17.3760 1.4480 0.4817 0.9320
## 5 19.9117 0.3378 0.5685 0.5385 0.2253 13.7559 2.7512 0.4358 0.5223
## Ptbiserial Frey McClain Dunn Hubert SDindex Dindex SDbw
## 2 0.8799 2.8318 0.2099 0.4676 0.0270 1.3933 0.5284 0.2281
## 3 0.6879 0.9602 0.4835 0.2546 0.0272 2.5155 0.4092 0.2077
## 4 0.6656 1.1769 0.5288 0.3593 0.0300 2.3116 0.3269 0.1564
## 5 0.6125 0.8156 0.6230 0.3593 0.0314 3.0478 0.2695 0.1515
##
## $All.CriticalValues
## CritValue_Duda CritValue_PseudoT2 Fvalue_Beale
## 2 0.1155 130.1288 0.2735
## 3 -0.1908 -37.4469 0.1507
## 4 -0.0027 -4016.9908 0.2566
## 5 -0.3258 -16.2785 0.1232
##
## $Best.nc
## KL CH Hartigan CCC Scott Marriot TrCovW TraceW
## Number_clusters 2.0000 2.0000 3.0000 5.0000 3.0000 3.000 4.0000 3.0000
## Value_Index 12.6385 105.1881 4.5336 2.5756 39.2084 38.954 11.9832 2.0891
## Friedman Rubin Cindex DB Silhouette Duda PseudoT2 Beale
## Number_clusters 3.0000 4.0000 5.0000 2.0000 2.0000 2.0000 2.000 2.0000
## Value_Index 15.4773 0.9625 0.3378 0.4436 0.7114 0.4121 24.252 1.3473
## Ratkowsky Ball PtBiserial Frey McClain Dunn Hubert
## Number_clusters 2.0000 3.0000 2.0000 2.0000 2.0000 2.0000 0
## Value_Index 0.6355 3.0553 0.8799 2.8318 0.2099 0.4676 0
## SDindex Dindex SDbw
## Number_clusters 2.0000 0 5.0000
## Value_Index 1.3933 0 0.1515
##
## $Best.partition
## 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016
## 1 1 1 2 2 2 2 2 2 2 2Tras utilizar diversos métodos para ver cuál es el número óptimo de grupos, la función concluye que el mejor número de clúster a formar es 2, de acuerdo a 13 índices. Por otra parte, 6 índices recomiendan que el número de grupos a formar sea 3. Vamos a ceñirnos a esta última conclusión y vamos a formar 3 clústeres, ya que si formamos 2 resultarían desproporcionados en observaciones.
4.3. Elección óptima de clústeres: Método no jerárquico.
Como hemos decidido formar el número de grupos, vamos a utilizar el método no jerárquico de k-medias. En primer lugar, fijamos la semilla para que el resultado no varíe, y posteriormente, creamos una variable que contenga 3 grupos diferenciados:
set.seed(1)
tresgrupos <- kmeans(scale(datos3), 3)y a continuación observamos el dendrograma y un gráfico de representación de los grupos:
factoextra::fviz_dend(cluster, 3)
factoextra::fviz_cluster(tresgrupos, datostip, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE)
En el dendrograma se observan diferentes colores en función del grupo al que pertenece cada observación. En este caso vemos como podamos más abajo. También lo podemos observar en el gráfico con mayor claridad, donde se ve un grupo claramente diferente (el verde), junto a dos grupos que se llegan a solapar en una parte (rojo y azul). También se observan los centroides de cada grupo.
Procedemos mediante k-medias a crear 3 grupos o clústeres:
kme <- kmeans(datostip, 3)
kme## K-means clustering with 3 clusters of sizes 5, 18, 4
##
## Cluster means:
## Renovables Emisiones_CO2
## 1 1.4599398 -1.7524689
## 2 -0.6314184 0.5970601
## 3 1.0164579 -0.4961845
##
## Clustering vector:
## 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016
## 2 2 3 3 3 3 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 0.6856463 5.2981353 1.2927588
## (between_SS / total_SS = 86.0 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Se han creado 3 grupos con 8, 13 y 6 observaciones. Más abajo se muestra a que grupo ha sido seleccionado cada observación.
datostipkm <- data.frame(datostip, GRUPO = factor(kme$cluster))
datos3$GRUPO <- factor(kme$cluster)
pander::pandoc.table(datos3)| Renovables | Emisiones_CO2 | GRUPO | |
|---|---|---|---|
| 1990 | 38.95 | 4545 | 2 |
| 1991 | 40.55 | 4487 | 2 |
| 1992 | 40.15 | 4347 | 2 |
| 1993 | 43.08 | 4267 | 2 |
| 1994 | 42.58 | 4250 | 2 |
| 1995 | 44.28 | 4301 | 2 |
| 1996 | 46.99 | 4409 | 2 |
| 1997 | 48.97 | 4323 | 2 |
| 1998 | 49.21 | 4322 | 2 |
| 1999 | 48.17 | 4265 | 2 |
| 2000 | 49.14 | 4294 | 2 |
| 2001 | 49.11 | 4364 | 2 |
| 2002 | 50.35 | 4337 | 2 |
| 2003 | 54.83 | 4433 | 2 |
| 2004 | 55.27 | 4450 | 2 |
| 2005 | 59.38 | 4434 | 2 |
| 2006 | 63.05 | 4446 | 2 |
| 2007 | 68.7 | 4404 | 2 |
| 2008 | 72.92 | 4304 | 3 |
| 2009 | 77.03 | 3951 | 3 |
| 2010 | 82.78 | 4071 | 3 |
| 2011 | 78.88 | 3929 | 3 |
| 2012 | 83.7 | 3868 | 1 |
| 2013 | 84.95 | 3781 | 1 |
| 2014 | 83.06 | 3614 | 1 |
| 2015 | 86.92 | 3652 | 1 |
| 2016 | 88.95 | 3637 | 1 |
4.4. Interpretación de las emisiones de CO2 y energías renovables en Europa.
Una vez tenemos 3 grupos relativamente diferentes, vamos a comentar que características tienen cada uno de ellos:
pander::pandoc.table(aggregate(cbind(Renovables, Emisiones_CO2) ~ GRUPO, data = datos3, FUN = mean))| GRUPO | Renovables | Emisiones_CO2 |
|---|---|---|
| 1 | 85.52 | 3710 |
| 2 | 49.6 | 4371 |
| 3 | 77.9 | 4064 |
El grupo 1 contiene los últimos 8 años (el último considerado es 2016), y tienen la característica de ser los años menos contaminantes y con mayor producción de energías renovables, lo cual es positivo de cara al futuro.
El grupo 2 contiene a las primeras observaciones (desde 1990 hasta 2002), y se caracterizan por ser los años con menor producción de renovables. Llama la atención que en estos años no se haya producido tanta contaminación como en el grupo intermedio, el cual es más avanzado en el tiempo.
El grupo 3 es el grupo intermedio entre los anteriores e incluye desde el 2003 hasta el 2008. Destaca que, a pesar del incremento de renovables, de media se contaminaron más estos años que en el grupo 2, a pesar de que en este último se produjeron menos renovables.
En los siguientes boxplot vemos al relación clara entre contaminación y producción de renovables (cuanto más renovables, menos contaminación y viceversa):
par(mfrow = c(1, 2))
boxplot(Renovables ~ GRUPO, main = "Boxplot de producción de renovables", col = c("palevioletred1","orange","green2"), data = datos3, las = 1)
boxplot(Emisiones_CO2 ~ GRUPO, main = "Boxplot de las emisiones de CO2", col = c("palevioletred1","orange","green2"), data = datos3, las = 1)
Es destacable que entre los años 2003-2008 las emisiones de CO2 aumentaron, a pesar que en el periodo anterior (1990-2002) la tendencia era a la baja. No es hasta 2009 cuando hay un salto muy claro en la bajada de emisiones de C02 (de 4304mtoe en 2008 a 3951mtoe en 2009).
También podría utilizarse el análisis de la varianza, o la MANOVA en el caso de análisis múltiple para comprobar la significación de los clústeres. En este análisis no se va a realizar para no repetir análisis similares a los del primer apartado.
Como conclusión, la tendencia actual en la producción de energías renovables es a la alza, al igual que es a la baja las emisiones de CO2, sin embargo, esto no es suficiente, ya que los expertos siguen afirmando que es demasiado poco el esfuerzo que se realiza en intentar evitar el cambio climático, la temperatura del planeta sigue subiendo y el permafrost se está descongelando, dando lugar a todavía más emisiones nocivas para el planeta. Debe tomarse conciencia y actuar lo más rápido posible para evitar catástrofes mayores.