Relación entre la agricultura y la temperatura mundial
Daniel Rodríguez García
7/5/2020
El cambio en el clima del planeta se debe a factores naturales y antropogénicos. Si bien es cierto, la temperatura ambiental fluctua con el tiempo, debido a factores externos e internos del planeta, el patrón se ha alterado a causa del hombre. Desde la revolución industrial, hubo una depredación desmedida de los recursos para sustentar la demanda de alimento y transporte del mundo. Por ello, se ha talado grandes áreas de bosque y se ha transformado en zonas industriales, agrícolas y urbanas.
Como se mencionó, las actividades agrícolas han sustituido áreas naturales de gran extensión. Esto puede haber provocado un incremento en la temperatura ambiental por diversos factores. En primer lugar, el uso de agroquímicos pueden generar gases de efecto invernadero que se acumulan en la atmósfera degradándola. También, las áreas cultivadas no tienen la misma capacidad de absorción de calor, como si lo tienen los bosques gracias a su composición estratigráfica. Por último, dependiendo del cultivo y la técnica utilizada, puede haber mayor reflectancia de los rayos solares, que asociado a los gases de efecto invernadero, pueden provocar un incremento de la energía calorífica dentro de la atmósfera. Esto provocaría variaciones en la temperatura ambiental.
En este trabajo se quiso determinar la relación que existe entre las áreas cultivadas y la variación en la temperatura. Para ello, se utilizaron técnicas de exploración de datos para ver el comportamiento de la variable Temperatura a lo largo de los años. También se eligieron los países con mayor área cosechada durante 1960 y el 2018 y los cultivos con mayor terreno abarcado para ver cómo variaban a lo largo del tiempo y cómo variaban entre sí. Por último se hizo una relación lineal para ver la asociación entre las áreas cosechadas y la variación de la temperatura. Esto nos permitió crear un modelo que predijera la variación en la temperatura dependiendo de las áreas cultivadas.
Muestreo aleatorio de los países
En primer lugar, se debe escoger una muestra de todos los países de forma aleatoria e independiente. Esto significa que se escogen países al azar de los 5 continentes y luego se obtienen datos de estos países.
Así se ve la base de datos. Podríamos decir que está “sucia” o con información no muy útil. Para eso, hacemos un arreglo y ordenamos la matriz para convertirla en algo más entendible.
| Continente | CONTINENTE | PAIS | CAPITAL |
|---|---|---|---|
| africa | Africa | ANGOLA | LUANDA |
| NA | Africa | ARGELIA | ARGEL |
| NA | Africa | BENIN | PORTO-NOVO |
| NA | Africa | BOTSUANA | GABERONES |
| NA | Africa | BURKINA FASO | UAGADUGU |
| NA | Africa | BURUNDI | BUYUMBURA |
Este arreglo selecciona las columnas que quiero, luego agrupo los datos por CONTINENTE para tenerlo más ordenado
Así queda después de realizar el arreglo y puedo ahora sí, hacer mi muestreo aleatorio.
| CONTINENTE | PAIS |
|---|---|
| Africa | ANGOLA |
| Africa | ARGELIA |
| Africa | BENIN |
| Africa | BOTSUANA |
| Africa | BURKINA FASO |
| Africa | BURUNDI |
Escogencia aleatoria por continente
Ya teniendo la base de datos lista, se hace un muestreo aleatorio de los países. Sin embargo, quiero tener 10 países por cada continente para hacer un muestreo aleatorio estratificado.
El set.seed(19) es una forma de decirle a R el lugar donde quiero empezar la aleatoriedad. Por lo que cada vez que haya un set.seed(19), los resultados serán los mismos.
Nota aclaratoria: Puedo poner el número que quiera, no es necesario solamente 19, pero el otro número va a arrojar otros resultados diferentes.
set.seed(19)
#Africa
Africa <- paises[1:55,] # Estas son las filas donde están los países africanos
muestraAfrica <- sample(1:55, 10, replace = F) # Muestreo sin reemplazo (Ninguno se repite)
Africa1 <- Africa[muestraAfrica,] # Selección de las filas del muestreo
#America
America <- paises[56:91,] # Filas de países americanos
muestraAmerica <- sample(1:35, 10, replace = F) # Son 35 países americanos, de los cuales ocupo solamente 10.
America1 <- America[muestraAmerica,]
#Asia
Asia <- paises[92:137,]
muestraAsia <- sample(1:45, 10, replace = F)
Asia1 <- Asia[muestraAsia,]
#Europa
Europa <- paises[138:184,]
muestraEuropa <- sample(1:50, 10, replace = F)
Europa1 <- Europa[muestraEuropa,]
#Oceania
Oceania <- paises[185:198,]
muestraOceania <- sample(1:14, 10, replace = F)
Oceania1 <- Oceania[muestraOceania,]Unificación de 5 data.frame diferentes
De esta forma, tenemos 5 bases de datos independientes para cada continente. Pero, la queremos unificada, entonces aplico:
Aquí tenemos los 50 países seleccionados.
## [1] "SEYCHELLES" "ZAMBIA" "BENIN"
## [4] "BURKINA FASO" "GUINEA" "Djibouti"
## [7] "SUAZILANDIA" "MOZAMBIQUE" "GABON"
## [10] "NIGER" "JAMAICA" "BARBADOS"
## [13] "HAITI" "CANADA" "NICARAGUA"
## [16] "ESTADOS UNIDOS" "EL SALVADOR" "BELICE"
## [19] "CUBA" "HONDURAS" "COREA DEL NORTE"
## [22] "MALASIA" "CHINA" "VIETNAM"
## [25] "ISRAEL" "JAPON" "TAIWAN"
## [28] "AFGANISTAN" "CATAR" "INDONESIA"
## [31] "MONACO" "LITUANIA" "AZERBAIYAN"
## [34] "ESLOVAQUIA" "NORUEGA" "UCRANIA"
## [37] "REINO UNIDO" "LETONIA" "PORTUGAL"
## [40] "ALBANIA" "PAPUA NUEVA GUINEA" "ISLAS MARSHALL"
## [43] "PALAOS" "VANUATU" "NAURU"
## [46] "FIYI" "KIRIBATI" "SAMOA"
## [49] "ISLAS SALOMON" "TUVALU"Exploración de los datos
Cultivos
Ahora, vamos a utilizar unos datos nuevos. Estos son datos del área cosechada de los cultivos reportados a la F.A.O. Sin embargo, solamente vamos a utilizar los 10 productos con más área cosechada.
La base de datos a utilizar tiene esta disposición:
| Continente | Ambito | Pais | Producto | Fecha | Hectareas |
|---|---|---|---|---|---|
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1961 | 600 |
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1962 | 600 |
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1963 | 600 |
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1964 | 600 |
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1965 | 600 |
| QC | Cultivos | Afganistan | Aceitunas, olivas | 1966 | 600 |
Comportamiento de los datos por año
Ahora, se puede hacer cálculos dependiendo de cómo agrupemos los datos. En este caso, se agruparon por Año (Fecha en la base de datos). Luego de agrupados se hace un cálculo del promedio de hectareas, el cuál se guarda en una columna llamada Area_media.
cultivos2 <- cultivos %>%
group_by(Fecha) %>%#Agrupo por la columna año
na.exclude(Fecha) %>% #Excluyo las filas con NA's (Si lo hiciera con na.omit, no me deja hacer el cálculo de la media)
summarise(Area_media = mean(Hectareas)) #Calculo el promedio de hectareas cosechadas por añoCuando vemos el comportamiento de la cantidad de área cosechada por año, hubo una disminución de este dato en los años 90’s.
library(RColorBrewer)
g1 <- ggplot(cultivos2, aes(x=Fecha, y = Area_media))+
geom_line()+
geom_point(aes(group = seq_along(Fecha)), size = 1)+
labs(x = "Año", y = "Area media cosechada")+
theme_classic()+
transition_reveal(Fecha)+
coord_cartesian(clip = 'off') +
enter_fade() +
exit_shrink()
animate(g1)
Comportamiento de los datos por año y por producto
Si se piensa en la diversificación de los productos por la demanda del mercado, no es raro suponer que mientras unos productos son poco atractivos con el tiempo, hay otros que aumentan su producción. Es por eso que se realizó una agrupación por fecha y por producto, luego se calculó el área promedio y se ordenó de forma que estuviera el producto y el año con mayor área cultivada de primero.
cultivos3 <- cultivos %>%
na.exclude() %>% #Excluyo las filas con NA's
group_by(Fecha,Producto) %>% #Agrupo por año y por producto
summarise(Area_media = mean(Hectareas))%>%
#Calculo el promedio de los datos por año y por cultivo (Los dos al mismo tiempo)
arrange(desc(Area_media)) # Ordeno de manera descendente los datos tomando como criterio el área media cosechadaAquí podemos ver los diez primeros registros:
| Fecha | Producto | Area_media |
|---|---|---|
| 1981 | Trigo | 4323903 |
| 1990 | Trigo | 4295966 |
| 1982 | Trigo | 4255138 |
| 1984 | Trigo | 4137671 |
| 1985 | Trigo | 4094501 |
| 1980 | Trigo | 4091214 |
| 1989 | Trigo | 4081612 |
| 1991 | Trigo | 4060725 |
| 1986 | Trigo | 4060104 |
| 1976 | Trigo | 4056529 |
Ahora, escojo los 10 primeros productos con mayor área cosechada. Esto se determina viendo cuáles son los primeros diez productos de la tabla anterior
cultivos3_seleccionados <- cultivos3 %>%
filter(Producto %in% c("Trigo", "Soja", "Maiz", "Colza", "Arroz, cascara", "Cebada", "Aceite, nuez de palma", "Sorgo","Caucho natural", "Mijo"))g2 <- ggplot(cultivos3_seleccionados,aes(x=Fecha,y =Area_media, col = Producto))+
geom_line(size = 1)+
geom_point(aes(group = seq_along(Fecha)), size = 1)+
scale_colour_brewer(palette="Set3")+
theme_classic()+
transition_reveal(Fecha)+
coord_cartesian(clip = 'off') +
enter_fade() +
exit_shrink()
animate(g2)
Cultivos que tienen similar comportamiento
La impresión que deja el gráfico anterior es que las áreas cosechadas de los diez productos más “importantes” no se comportan de manera similiar a lo largo del tiempo. Para saber cuáles si son similares entre sí, podemos hacer un análisis de componentes principales (PCA). Este análisis disminuye las dimensiones de los datos. Lo cual quiere decir que en un gráfico donde hay muchas variables (En este caso productos) y por ende muchas dimensiones. Se puede acomodar para verlo en 2D, como la pantalla de nuestra computadora. Este tipo de análisis necesita una configuración de los datos diferente:
## Warning: Missing column names filled in: 'X1' [1]| X1 | Aceite, nuez de palma | Arroz,cascara | Caucho natural | Cebada | Colza | Maiz | Mijo | Soja | Sorgo | Trigo |
|---|---|---|---|---|---|---|---|---|---|---|
| 1961 | 383993.5 | 1670824 | 582993.0 | 1121001 | 244966.4 | 1417851 | 990735.8 | 2242330 | 885096.8 | 3607792 |
| 1962 | 359356.1 | 1724133 | 584240.0 | 1225118 | 213019.5 | 1338390 | 922534.6 | 2218091 | 904078.4 | 3388246 |
| 1963 | 372837.6 | 1714216 | 605420.0 | 1167708 | 224463.5 | 1440915 | 995495.5 | 2262811 | 944537.3 | 3419183 |
| 1964 | 373999.8 | 1810895 | 603556.6 | 1068820 | 282331.1 | 1439794 | 955417.9 | 2388540 | 918697.9 | 3668605 |
| 1965 | 387268.0 | 1827362 | 637560.0 | 1047970 | 315106.0 | 1414964 | 958783.8 | 2398347 | 961175.2 | 3591346 |
| 1966 | 376453.6 | 1864373 | 615150.2 | 1093653 | 307339.9 | 1507755 | 947368.0 | 2464532 | 910262.7 | 3564645 |
Hacemos unos pequeños arreglos a la base de datos para poder meterlo en el análisis:
filas <- cultivos_pca$X1 #Esta es una columna donde está el nombre de los años
cultivos_pca1 <- cultivos_pca[,-1] #Quito dicha columna de la base
cultivos_pca1 <- as.matrix(cultivos_pca1) #transformo a matriz
rownames(cultivos_pca1) <- filas #Nombro las filas con el año Este cambio permitirá ver en el gráfico los años como puntos y los productos como líneas
## Aceite, nuez de palma Arroz,cascara Caucho natural Cebada Colza
## 1961 383993.5 1670824 582993.0 1121001 244966.4
## 1962 359356.1 1724133 584240.0 1225118 213019.5
## 1963 372837.6 1714216 605420.0 1167708 224463.5
## 1964 373999.8 1810895 603556.6 1068820 282331.1
## 1965 387268.0 1827362 637560.0 1047970 315106.0
## 1966 376453.6 1864373 615150.2 1093653 307339.9
## Maiz Mijo Soja Sorgo Trigo
## 1961 1417851 990735.8 2242330 885096.8 3607792
## 1962 1338390 922534.6 2218091 904078.4 3388246
## 1963 1440915 995495.5 2262811 944537.3 3419183
## 1964 1439794 955417.9 2388540 918697.9 3668605
## 1965 1414964 958783.8 2398347 961175.2 3591346
## 1966 1507755 947368.0 2464532 910262.7 3564645Para hacer el PCA ocupamos instalar y leer estos paquetes:
library(devtools)
library(factoextra) # este se tiene que instalar con el siguiente código devtools::install_github("kassambara/factoextra")Utilizamos la siguiente función con el argumento scale = T que homogeniza la varianza en 1.
Luego de hecho el análisis, lo vamos a graficar de dos diferentes maneras. El primero es nada más la disposición de las variables entre sí. Están coloreadas dependiendo al grado de contribución para hacer los componentes principales.
fviz_pca_var(cult_pca,
col.var = "contrib", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, # Avoid text overlapping
title = "Disposición de las variables"
)
Ahora podemos ver cómo se distribuyen los países entre los productos. Podemos ver países que se acercan a productos como el Maíz, el Sorgo y la Cebada. Esto puede inferir que en esos países, el cultivo más importante es al cuál más se acercan.
fviz_pca_biplot(cult_pca, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969", # Individuals color
alpha.ind = 0.5,
title = "Biplot entre años y productos"
)
Áreas cosechadas por países
Esta base de datos es la que tiene solamente los datos de los países seleccionados aleatoriamente.
En primer lugar, debemos extraer la información necesaria. Se seleccionan las columnas que necesito, filtro por los productos principales y se agrupa por país. Luego se calcula el área media cosechada y se ordena de forma descendente.
cultivos4 <- cultivos %>%
select(Fecha, Pais, Hectareas, Producto) %>%
filter(Producto %in% c("Trigo", "Soja", "Maiz", "Colza", "Arroz, cascara", "Cebada", "Aceite, nuez de palma", "Sorgo","Caucho natural", "Mijo")) %>%
na.exclude() %>%
group_by(Pais) %>%
summarise(Area_media = mean(Hectareas))%>%
arrange(desc(Area_media))# Ordena la base con los primeros cultivos en importancia| Pais | Area_media |
|---|---|
| Estados Unidos de América | 11095715.7 |
| China, Continental | 11075973.8 |
| Canada | 4013935.8 |
| Indonesia | 3744964.0 |
| Nigeria | 2106287.8 |
| Ucrania | 1694111.3 |
| Viet Nam | 1649389.3 |
| Reino Unido | 1116148.0 |
| Malasia | 953799.5 |
| Burkina Faso | 589896.2 |
Ahora, se filtra la base por los primeros 10 países con mayor área cosechada promedio.
cultivos4_seleccionados <- cultivos %>%
filter(Pais %in% c("Estados Unidos de América","China, Continental", "Canada","Indonesia","Nigeria","Ucrania", "Viet Nam", "Reino Unido", "Malasia","Burkina Faso")) %>%
select(Fecha, Pais, Hectareas, Producto) %>%
filter(Producto %in% c("Trigo", "Soja", "Maiz", "Colza", "Arroz, cascara", "Cebada", "Aceite, nuez de palma", "Sorgo","Caucho natural", "Mijo")) %>%
na.exclude() %>%
group_by(Fecha,Pais,Producto) %>%
summarise(Area_media = mean(Hectareas))Se puede visualizar el comportamiento del dato área media cosechada por país, pero agrupado por producto.
g3 <- ggplot(cultivos4_seleccionados,aes(x = Fecha, y =Area_media,col = Pais))+
geom_line(size = 1)+
geom_point(aes(group = seq_along(Fecha)), size = 1)+
theme(axis.text.x = element_text(angle = 45))+
theme_classic()+
scale_colour_brewer(palette="Set3")+
facet_wrap(~Producto,scales = "free_y")+
theme(legend.position = "bottom")+
transition_reveal(Fecha)+
coord_cartesian(clip = 'off') +
enter_fade() +
exit_shrink()
animate(g3)
Esto nos da una idea del comportamiento de las cosechas por país a lo largo del tiempo.
Variación de la Temperatura
El calentamiento global y el cambio climático se debe en parte a las acciones del ser humano. Tanto las industrias como las actividades agropecuarias, han provocado un desbalance en la cantidad de gases de efecto invernadero y un aumento en la temperatura. En primer lugar podemos ver cómo ha variado la temperatura a lo largo de los años, con respecto a una temperatura de referencia (temperatura promedio entre los años 40 y 60).
Se filtran los datos por los países con más área media cosechada.
temperatura1 <- temperatura %>%
select(Variacion_temperatura,Fecha,Pais) %>%
filter(Fecha != 2019) %>%
na.exclude() %>%
filter(Pais %in% c("Estados Unidos de América","China, Continental", "Canada","Indonesia","Nigeria","Ucrania", "Viet Nam", "Reino Unido", "Malasia","Burkina Faso"))Vemos que la variación de la temperatura aumenta gradualmente con respecto al tiempo. Esto demuestra que hay un factor que está alterando la temperatura en los países escogidos.
g4 <- ggplot(temperatura1,aes(x = Fecha, y =Variacion_temperatura,col = Pais))+
geom_line(size = 1)+
geom_point(aes(group = seq_along(Fecha)), size = 1)+
theme(axis.text.x = element_text(angle = 45))+
theme_classic()+
scale_colour_brewer(palette="Set3")+
labs(x= "Año", y="Variación de la temperatura (° C)")+
theme(legend.position = "bottom")+
transition_reveal(Fecha)+
coord_cartesian(clip = 'off') +
enter_fade() +
exit_shrink()
animate(g4)
Relación entre área cosechada y la variación de la temperatura
Ahora, para hacer una relación entre estas dos variables, hay que tenerlas en un sólo set de datos. Esto no lo podemos hacer como se unieron los paises al principio porque no tienen las mismas columnas. Por ello se debe usar otra estrategia:
cultivos5 <- cultivos %>%
filter(Pais %in% c("Estados Unidos de América","China, Continental", "Canada","Indonesia","Nigeria","Ucrania", "Viet Nam", "Reino Unido", "Malasia","Burkina Faso")) %>%
select(Fecha, Pais, Hectareas, Producto) %>%
na.exclude() %>%
group_by(Pais,Fecha) %>%
summarise(Area_media = mean(Hectareas))Creamos una columna nueva en cultivos5 llamada temperatura que sea igual a la columna variación de la temperatura de la base de datos temperatua1.
| Pais | Fecha | Area_media | temperatura |
|---|---|---|---|
| Burkina Faso | 1961 | 110578.1 | -0.302 |
| Burkina Faso | 1962 | 118392.6 | -0.251 |
| Burkina Faso | 1963 | 121130.3 | 0.126 |
| Burkina Faso | 1964 | 124244.1 | -0.231 |
| Burkina Faso | 1965 | 124768.4 | -0.135 |
| Burkina Faso | 1966 | 128527.8 | -0.149 |
Luego, hacemos el cálculo del promedio por año del área cultivada y de la temperatura. Así, tenemos el comportamiento entre las dos variables con datos de promedio anual.
cultivos_fecha_temperatura <- cultivos5 %>%
group_by(Fecha) %>%
summarise(Area_media = mean(Area_media),
Temperatura = mean(temperatura))Ahora, para ver la asociación que hay entre el área cosechada y la variación de la temperatura debemos utilizar la estadística. En este caso, por ser dos variables numéricas, utilizamos una regresión lineal.
ggplot(cultivos_fecha_temperatura,aes(x=Area_media,y=Temperatura))+
geom_line()+
geom_smooth(method = "lm", se = F)+
theme_classic()+
labs(x = "Area media cosechada", y= "Variación de la temperatura (°C)" )## `geom_smooth()` using formula 'y ~ x'
Observamos un aumento de la variación de la temperatura mientras aumenta el área media cosechada (R2 = 0.29, Estimate = 8.61 * 10-6, p-value < 0.005).
modelo_lineal1 <- lm(Temperatura~Area_media,data = cultivos_fecha_temperatura)
summary(modelo_lineal1)##
## Call:
## lm(formula = Temperatura ~ Area_media, data = cultivos_fecha_temperatura)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70281 -0.32084 0.00567 0.33960 0.77605
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.267e+00 9.773e-01 -4.366 5.51e-05 ***
## Area_media 8.619e-06 1.780e-06 4.843 1.05e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4015 on 56 degrees of freedom
## Multiple R-squared: 0.2952, Adjusted R-squared: 0.2826
## F-statistic: 23.45 on 1 and 56 DF, p-value: 1.05e-05Predicciones de la variación de la temperatura dependiendo del área
Con estos datos se puede predecir la variación de la temperatura en años siguientes, solamente sabiendo el área media cultivada. Para eso, utilizamos los modelos lineales generalizados.
Para hacer estas predicciones utilizamos el modelo que ya construimos como apoyo. Colocamos las nuevas predicciones como una nueva columna.
cultivos_fecha_temperatura$prediccion <- predict(modelo_lineal2, type = "response") #Lo coloco en una columna dentro del data.frame
head(cultivos_fecha_temperatura)| Fecha | Area_media | Temperatura | prediccion |
|---|---|---|---|
| 1961 | 510558.6 | 0.0245556 | 0.1335256 |
| 1962 | 512510.7 | -0.2092222 | 0.1503503 |
| 1963 | 516243.5 | -0.0044444 | 0.1825233 |
| 1964 | 527558.2 | -0.1142222 | 0.2800433 |
| 1965 | 531498.7 | -0.2706667 | 0.3140062 |
| 1966 | 534268.8 | 0.0164444 | 0.3378818 |
En este caso, la línea roja viene siendo las predicciones hechas por el modelo y la negra los datos reales (obtenidos en el muestreo).
g6 <- ggplot(cultivos_fecha_temperatura,aes(x=Fecha,y=Temperatura))+
geom_line()+
geom_line(aes(x = Fecha, y = prediccion),color= "Red" )+
geom_point(aes(group = seq_along(Fecha)), size = 1)+
theme_classic()+
transition_reveal(Fecha)+
coord_cartesian(clip = 'off') +
enter_fade() +
exit_shrink()
animate(g6)
Sin embargo, este modelo no se ajusta bien a los datos, porque el R2 muestra que solamente el 30% se asocian en tre sí. Además, el valor de Gini es lejano a 1 (0.44). Esto me dice que el modelo lineal simple no nos sirve para predecir, ocupo otras estrategias de Machine Learning (Random Forest, Árboles de decisión,… ).
library(WVPlots)
GainCurvePlot(cultivos_fecha_temperatura,"prediccion","Temperatura", title= "Ajuste del modelo")