Aplicación de modelos mixtos para ajustar el volumen de producción de cultivos de café en Antioquia

# Paquete utilizado para el ajuste de modelos mixtos.
library(nlme)

# Paquete utilizado para realizar análisis descriptivo.
library(tidyverse)

# Paquete para el análisis de residuales
library(predictmeans)

# Paquete para gráficos
library(ggplot2)
library(ggpubr)

1. Descripción

La agricultura es una de las principales actividades económicas del departamento de Antioquia y el café es uno de los productos mas cultivados. Son 80.238 productores, de los 550 mil que tiene el país, que todos los días se levantan con sus hijos, esposas y trabajadores a recorrer 104 mil fincas que se extienden en 2.685 veredas de 93 de los 125 municipios del departamento. Adicionalmente es uno de los productos que representan un alto porcentaje de participación en la exportación, teniendo esto en cuenta en importante implementar modelos que expliquen el comportamiento de la producción de café en el departamento teniendo en cuenta variables como el área de producción y el municipio de cultivo.

1.1 Lectura de la base de datos

Los datos fueron tomados de la página web Datos Abiertos Colombia, el link de descarga es https://www.datos.gov.co/Agricultura-y-Desarrollo-Rural/Areas-cultivadas-y-produccion-agr-cola-en-Antioqui/t2ca-uae5https://www.datos.gov.co/Agricultura-y-Desarrollo-Rural/Areas-cultivadas-y-produccion-agr-cola-en-Antioqui/t2ca-uae5

agricola <- read.csv("Agricola.csv", encoding = "UTF-8")

El conjunto de datos contiene información de los cultivos en Antioquia desde el año 1990, entendiendo que la agricultura ha evolucionado se decide seleccionar para el análisis y ajuste de modelos la información desde el año 2014 hasta 2020 y únicamente los cultivos de café ya que es uno de los productos que se exporta en gran cantidad desde el departamento.

Selección de columnas relevantes, año y cultivo a analizar.

data <- agricola %>% 
  filter(Año %in% c("2014", "2015", "2016", "2017", "2018", "2019", "2020"))
cafe <- data %>% filter(str_detect(Rubro, "Café")) %>% 
  select(c("Municipio", "Área.Producción", "Volumen.Producción"))

1.2 Descripción de variables

La base de datos contiene 93 municipios del departamento de Antoquia donde se cultiva café y 634 registros uno por cada combinación de municipio y área de producción y se observa el volumen en toneladas producido.

1.3 Vista del conjunto de datos

2. Objetivos

1. Encontrar el mejor modelo predictivo que explique la producción en toneladas de los cultivos de Café en el departamento de Antioquia entre los años 2014 a 2020 discriminando por los diferentes municipios del departamento, además de observar la significancia de los efectos fijos y aleatorios mediante pruebas de hipótesis para finalmente escoger el mejor modelo.

2. Exhibir la utilidad e importancia de los modelos mixtos cuando se tienen datos agrupados.

3. Análisis descriptivo

Antes de iniciar con el ajuste de modelos es importante realizar un análisis descriptivo que nos permita visualizar las relaciones existentes entre la variable objetivo y las variables predictoras, en este caso se quiere observar el comportamiento del volumen de producción, para esto se crea un gráfico de dispersión.

ggplot(data = cafe, aes(x = Área.Producción, y = Volumen.Producción, 
                           color = Municipio)) +
  geom_point() +
  theme_bw() + 
  labs(y = "Volumen de Producción") +
  ggtitle ("Gráfico de dispersión de Área de Producción vs Volumen de Producción")

Se observa en la nube de puntos que el volumen de producción de los municipios de Antioquia crece conjuntamente cada vez que el área de producción se hace más grande. Además, parece que la variabilidad aumenta a medida que las áreas de producción son mas grandes, lo cual indica que ajustar un modelo de regresión lineal clásico presentaría problemas con el supuesto de varianza constante.

A continuación, se presenta un gráfico de dispersión por cada municipio, para analizar de una mejor manera el comportamiento de la producción en cada municipio.

ggplot(data = cafe, aes(x = Área.Producción, y = Volumen.Producción, 
                           color = Municipio)) +
  geom_point() +
  theme_bw() + 
  labs(y = "Volumen de Producción", x = "Área Producción") +
  ggtitle ("Gráfico de dispersión de Área de Producción vs Volumen de Producción") +
  facet_wrap(~ Municipio) +
  theme(legend.position = "none") 

En la figura anterior se observa que la nube de puntos es creciente para todos los municipios en los que se cultiva café, existen algunos municipios como Andes, Betania, Betulia, Cuidad Bolivar, Concordia y Salgar que poseen valores muy altos tanto de áreas cultivadas como de producción, esto nos indica que estos son los municipios que posiblemente mayores exportaciones de café generan.

También se observa que los interceptos son diferentes en cada municipio, por lo que la creación de modelos mixtos puede ser beneficiosa para explicar el comportamiento del volumen de producción en función del área cultivada.

4. Modelos a considerar

1. Modelo de regresión lineal clásico.

\[produccion_j \sim N(\mu, \sigma^2)\]

\[\mu_j = \beta_0 + \beta_1AreaProduccion_j\]

2. Modelo de regresión lineal mixto con intercepto aleatorio.

\[produccion_{ij} \sim N(\mu_{ij}, \sigma^2)\] \[\mu_{ij} = \beta_0 + \beta_1 AreaProduccion_{ij} + \beta_{0i}\]

\[b_{0} \sim N(0, \sigma^2_{b0})\]

\[i = 1,~2,...,~93\]

3. Modelo de regresión lineal mixto con interccepto y pendiente aleatoria.

\[produccion_{ij} \sim N(\mu_{ij}, \sigma^2)\]

\[\mu_{ij} = \beta_0 + \beta_1 AreaProduccion_{ij} + \beta_{0i} + \beta_{1i} ~AreaProduccion_{ij}\]

\[ {b_0 \choose b_1} \sim N \left({0 \choose 0},{\sigma^2_{b0} \ \ \sigma^2_{b01} \choose \sigma^2_{b01} \ \ \sigma^2_{b1}} \right)\]

\[i = 1,~2,...,~93\]

5. Ajuste de modelos

A continuación se presenta el código utilizado para el ajuste de los modelos propuestos en la sección 4, para esto es utilizado el paquete \(nmle\) que esta orientado a la creación de modelos lineales mixtos con respuesta distribuida normal.

# Modelo 1
fit1 <- gls(Volumen.Producción ~ 1 + Área.Producción, data = cafe)

# Modelo 2
fit2 <- lme(Volumen.Producción ~ 1 + Área.Producción, 
            random =~ 1 | Municipio,
             data = cafe)

# Modelo 3
fit3 <- lme(Volumen.Producción ~ 1 + Área.Producción, 
            random=~ 1 + Área.Producción | Municipio,
             data = cafe)

6. Comparación de modelos

6.1 Prueba de hipótesis sobre componente de varianza del intercepto aleatorio

Se quiere determinar si introducir una variable de agrupación en el intercepto del modelo es significativa, para esto la hipótesis a probar es

\[H_0 : \sigma^{2}_{b0} = 0 ~~~~~~ vs ~~~~~~ H_A:\sigma^{2}_{b0} > 0\] con estadístico de prueba

\[lrt=−2×(logLik(modelo 1)−logLik(modelo 2))\] donde \(lrt \sim \chi_1^2\), se rechaza \(H_0\) con un nivel de significancia \(\alpha = 0.05\) si \(P(\chi_1^2 >lrt) < \alpha\). Para realizar esta prueba de hipótesis sobre esta componente de varianza es posible utilizar una tabla anova, sin embargo esta es una prueba conservativa, por lo que se debe tener cuidado con \(Valores~P\) grandes.

anova1 <- anova(fit1, fit2)
anova1 %>% mutate(call = NULL)

De la tabla anterior se tiene que \(lrt=421.4642\) y \(P(\chi_1^2 >lrt)<0.0001\) por lo que se rechaza la hipótesis nula y se tiene que la componente de varianza para el intercepto aleatorio es significativa y contribuye a la explicación del volumen de producción de los cultivos de café en Antioquia.

6.2 Prueba de hipótesis sobre componente de varianza de la pendiente aleatoria y componente de covarianza entre pendientes e interceptos aleatorios

Luego de identificar que la inclusión de interceptos aleatorios es significativo para el modelo es momento de realizar la prueba sobre la componente de varianza de las pendientes aleatorias y covarianza entre pendientes e interceptos aleatorios. La hipótesis a probar es

\[H_{0}: {\sigma^2_{b0} \ \ \sigma^2_{b01} \choose \sigma^2_{b01} \ \ \sigma^2_{b1}} = {\sigma^2_{b0} \ \ 0 \choose 0 \ \ ~~0} ~~~ vs ~~~ H_{A}: {\sigma^2_{b0} \ \ \sigma^2_{b01} \choose \sigma^2_{b01} \ \ \sigma^2_{b1}} = {\sigma^2_{b0}~~~~~~ \ \ \sigma^2_{b01} \neq 0 \choose \sigma^2_{b01} \neq0 \ \ \sigma^2_{b1}>0 }\] con estadístico de prueba

\[lrt= -2 \times (logLik(mododelo ~2) - logLik(modelo~3))\]

donde \(lrt \sim \chi^2_2\), se rechaza \(H_0\) con un nivel de significancia \(\alpha = 0.05\) si \(P(\chi_2^2 >lrt) < \alpha\). Para realizar esta prueba de hipótesis sobre esta componente de varianza y covarianza es posible utilizar una tabla anova, sin embargo esta es una prueba conservativa, por lo que se debe tener cuidado con \(Valores~P\) grandes.

anova2 <- anova(fit2, fit3)
anova2 %>% mutate(call = NULL)

De la tabla anterior se tiene que \(lrt=156.7391\) y \(P(\chi_2^2>lrt)<0.0001\) por lo que se rechaza la hipótesis nula y se tiene que la componente de varianza para la pendiente aleatoria, además de la componente de covarianza entre el intercepto y la pendiente aleatoria son significativas, por lo tanto contribuyen a la explicación de las toneladas de café producidas en el departamento de Antioquia.

6.3 Correlación entre valores ajustados y valores observados

correlaciones <- data.frame(cor_fit1 = cor(fitted(fit1),cafe$Volumen.Producción),
                            cor_fit2 = cor(fitted(fit2),cafe$Volumen.Producción),
                            cor_fit3 = cor(fitted(fit3),cafe$Volumen.Producción))
correlaciones

En la tabla anterior se presenta la correlación entre los valores predichos y reales para los tres modelos considerados, se observa que todos obtuvieron un buen desempeño, sin emargo, la mayor correlación la presenta el modelo 3, el cual tiene pendiente e intercepto aleatorio.

6.4 Criterios AIC y BIC

anova3 <- anova(fit1, fit2, fit3)
anova3 %>% mutate(call = NULL)

La tabla anterior indica que el mejor desempeño en cuanto a las medidas de AIC y BIC lo obtiene el modelo 3, ya que los valores de estas medidas son menores para este modelo.

Finalmente, de todas las tablas presentadas anteriormente y los criterios considerados es posible concluir que el mejor modelo es el número 3, el cual incluye pendientes e interceptos aleatorios, se procede a realizar el análisis de residuales para este modelo.

7. Análisis de residuales

7.1 Gráficos de residuales

En el gráfico de Residuales vs Valores ajustados, se observa un patrón similar a un embudo lo que nos indica una posible carencia de ajuste del modelo y posible violación del supuesto de varianza constante ya que la varianza aumenta hacia el final del eje x.

El gráfico Normal QQPlot para los residuales permite evidenciar que el supuesto de normalidad al parecer tiene problemas, ya que en los extremos varios puntos se alejan de la recta de normalidad, lo mismo se puede evidenciar con el histograma, los residuales estan centrados alrededor del cero pero tiene colas pesadas a ambos extremos.

7.2 Pruebas para los efectos aleatorios

Para la normalidad de los efectos \(b0\) y \(b1\) se puede observar que hay problemas en el supuesto de que estos distribuyen normal, pues hay varios puntos que están alejados de la linea (para ambos efectos), se observan colas pesadas.

7.3 Modelo con respuesta transformada

La manera más sencilla de solucionar la carencia de ajuste y varianza no constante evidenciada en los gráficos de residuales es aplicar una transformación logarítmica a la variable respuesta, es decir, a la producción en toneladas, ajustar de nuevo el modelo con el volumen de producción en escala logarítmica y posteriormente aplicarle la transformación inversa para obterner los valores ajustados.

7.3.1 Ajuste del modelo

# Aplicar logaritmo a la variable respuesta.
cafe$log_produccion <- log(cafe$Volumen.Producción)

# Modelo 3 con el volumen de produccion transformado
new_model <- lme(log_produccion ~ 1 + Área.Producción, 
            random=~ 1 + Área.Producción | Municipio,
             data = cafe)

# Agregar valores ajustados a la base de datos
cafe$pred_new <- exp(predict(new_model, type = 'response', newdata = cafe))

7.3.2 Gráficos de residuales

En el gráfico de Residuales vs Valores ajustados, ya no se observan patrones obvios, lo que nos indica que la transformación aplicada mejora notablemente el modelo, tanto en ajuste como en varianza constante.

El gráfico Normal QQPlot para los residuales permite evidenciar que el supuesto de normalidad presenta problemas, sin embargo las colas ya no son tan pesadas como se envidenciaba en el modelo 3.

7.3.3 Pruebas para los efectos aleatorios

Los supuestos de normalidad para los efectos \(b0\) y \(b1\) siguen con problemas, pero en general la transformación aplicada mejoro el ajuste del modelo.

8. Mejor modelo

Mediante las pruebas de hipótesis se concluye que el mejor modelo para ajustar la producción de café en Antioquia es el modelo 3. Pero en el análisis de residuales se encontró un problema con el supuesto de normalidad sobre la el volumen de producción. Por lo tanto, se usa una trasnformación logarítmica sobre el volumen de producción y se crea un modelo sobre esta variable transformada para corregir los supuestos inválidos del modelo 3.

Finalmente el mejor modelo ajustado es el siguiente:

\[Log[~ \widehat{\mu}_{ij}~]= 5.666495 + 0.001775~ ÁreaProducción_{ij} + \tilde{b}_{0i} + \tilde{b}_{1i} ~ÁreaProducción_{ij}\] \[con ~~\widehat{\mu}_{ij} = \widehat{E}[Volumen~de~Producción_{ij}~]\]

\[ {b_0 \choose b_1} \sim aprox ~ N \left({0 \choose 0},{1.803 \ \ \ -0.0021 \choose -0.0021 \ \ 2.567\times10^{-6}} \right)\]

\[i = 1,~2,...,~93\]

8.1 Representación gráfica de los modelos ajustados

Mejor modelo

Modelo 1

Modelo 2

Modelo 3

Referencias

Freddy Hernández Barajas, Jorge Leonardo López Martínez (2020-03-15). Modelos Mixtos con R. https://fhernanb.github.io/libro_modelos_mixtos/

Andrew Gelman, Jennifer Hill (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models.

Cómite de cafeteros de Antioquia. https://antioquia.federaciondecafeteros.org/cafe-de-antioquia/