Uno de los propósitos a nivel mundial es el de estimar la captura de dióxido de carbono uno principales gases causantes del calentamiento global y por tanto del cambio climático. Un grupo de investigadores están interesados en poder construir modelos que permitan la valoración de estos beneficio a través de información recogida sobre características de los árboles de una región húmeda, que les permita una estimación de la biomasa y así facilitar la toma de decisiones y la generación de políticas públicas. Como no es posible obtener el valor del peso del árbol sin cortar, se plantea la opción de estimarla a través de variables como la altura y el diámetro del tronco, información requerida para la estimación del valor de la biomasa.

\(*\) Se requiere ayude a los investigadores en su propósito utilizando la información contenida en la base de datos arboles suministrada.

\(*\) Proponga un modelo de regresión lineal simple que permita predecir el peso del árbol en función de las covariables que considere importantes y seleccionándolas de acuerdo con un proceso adecuado.

\(*\) Tenga en cuenta realizar una evaluación de la significancia de los parámetros, validación de los supuesto e interpretación de los resultados.

\(*\) Proponga un método de evaluación para los modelos estimados.

Carga de Datos

data("arboles")
head(arboles)
## # A tibble: 6 × 4
##      id  peso diametro altura
##   <int> <dbl>    <dbl>  <dbl>
## 1     1 13.7       4.7    5  
## 2     2 14.6       5.3    5.6
## 3     3 15.9       4.8    5.8
## 4     4  8.99      3.2    4.3
## 5     5  6.99      2.2    3.3
## 6     6 19.3       6.3    7.9

Descripción de las variables

summarytools::descr(arboles[,2:4])
## Descriptive Statistics  
## arboles  
## N: 90  
## 
##                     altura   diametro     peso
## ----------------- -------- ---------- --------
##              Mean     6.63       5.45    18.77
##           Std.Dev     1.80       1.45     8.16
##               Min     3.30       2.20     5.98
##                Q1     5.20       4.50    13.61
##            Median     6.45       5.40    17.48
##                Q3     7.90       6.50    22.82
##               Max    11.30       8.80    47.87
##               MAD     1.93       1.48     6.35
##               IQR     2.65       1.97     9.16
##                CV     0.27       0.27     0.43
##          Skewness     0.38      -0.07     1.11
##       SE.Skewness     0.25       0.25     0.25
##          Kurtosis    -0.44      -0.36     1.55
##           N.Valid    90.00      90.00    90.00
##         Pct.Valid   100.00     100.00   100.00

En la Base de Datos de los “Árboles”, encontramos a Y que es la variable de respuesta representada por el “Peso”, y X, la variable predictora que, dependiendo del caso, puede ser “Diámetro o Altura”.
Para el caso, es necesario recordar algunas definiciones que permitirá la fácil lectura en comparativo de las tres variables:

\(*\) Media es el promedio de todos los valores.

\(*\) Mediana es el valor central.

\(*\) La desviación estándar se mide con respecto a cero, la cual muestra la dispersión de la distribución de los datos.

\(*\) El Coeficiente de Variación (CV), es una medida de dispersión que permite el análisis de las desviaciones de los datos con respecto a la media y, de igual forma, las dispersiones que tienen los datos dispersos entre sí.

\(*\) La desviación media absoluta (MAD), es la distancia promedio entre cada valor y el promedio.

\(*\)La desviación media absoluta es una manera de describir la variación en un conjunto de datos.

\(*\) IQR es el rango intercuartílico.

\(*\) El skewness mide la asimetría de una distribución en relación con su media.

\(*\) El sesgo puede tomar valores negativos y positivos. Uno negativo indica que la cola está en el lado izquierdo de la distribución, que se extiende hacia valores más negativos. Un sesgo positivo indica que la cola está en el lado derecho de la distribución, que se extiende hacia valores más positivos.

\(*\) La kurtosis indica la altura y filo del pico central en la distribución de los datos, a mayor sea éste, más alto y puntiagudo será el pico central de los datos, esto es dado de acuerdo a un histograma.

Desarrollo del ejercicio

#G=ggarrange( g1, g2, ncol = 1, nrow = 2)
#g1=ggplot(arboles, aes(x=peso))+geom_density(fill="thistle1", color="gray8", alpha=0.7)

data <- data.frame(arboles$id,arboles$diametro,arboles$altura,arboles$peso) 
data_long <- melt(data, id = "arboles.id") 
cols <- c("violetred", "violet","thistle1")

ggplot(data_long, aes(x = value, colour = variable, fill = variable))+geom_density(alpha = 0.7,color="gray8")+scale_fill_manual(values = cols)+theme(legend.position="bottom")

options(repr.plot.width = 20, repr.plot.height =6)

\(*\) Posterior a la representación univariada de la base de datos, realizaremos un análisis de Regresión lineal, la cual, implica más que ajustar una linea recta a las observaciones de las variables, se podria pensar que, se requiere de tres pasos principales para realizar una buena estimación.

\(*\) La primera consiste en conocer el comportamiento las varibales (diámetro y altura) con respecto a la variable dependiente, en este caso es el peso de los árboles, es decir, la correlacion y dirección de los datos, de manera que se pueda seleccionar las mejores variables independientes para realizar la estimación.

\(*\) Posteriormente, se realiza la estimación del modelo.

\(*\) Luego, se requiere validar el modelo (ver si los supuestos se cumplen), es fundamental que esto se dé, para que el análisis e inferencia que se realice a partir de los resultados del modelo sea correcta (insesgada, consistete, y eficiente).

\(*\) Dado con lo anterior, es posible realizar análisis estadístico y pronóstico.

Estos son los grandes pasos que se van a evidenciar en este trabajo. Inicialmente empezaremos con el estudio de las correlaciones de las dos variables, altura y diámetro con respecto al peso de los árboles. 

scatterplot(arboles$diametro, arboles$peso,            # Datos
            pch = 19,        # Símbolo pch
            col = 1,         # Color de los puntos
            smooth = FALSE,  # Sin regresión suavizada
            regLine = FALSE, main=c("Correlacion",cor.test(arboles$diametro, arboles$peso)$estimate)) # Sin la regresión lineal

scatterplot(arboles$altura, arboles$peso,            # Datos
            pch = 19,        # Símbolo pch
            col = 1,         # Color de los puntos
            smooth = FALSE,  # Sin regresión suavizada
            regLine = FALSE, main=c("Correlacion",cor.test(arboles$altura, arboles$peso)$estimate)) # Sin la regresión lineal

Con estos resultados podemos ver que ambas variables tienen una correlacion alta y positiva con la variable dependiente (peso), sin embargo, se puede evidenciar que la variable diámetro, posee una mayor correlacion (0.91), que la altura (0.85), por lo tanto, tomaremos el diámetro como variable independiente (explicativa).

Iniciaremos estimando un modelo lineal sin transformaciones, sin embargo, en las correlaciones de las variables, asi como en la distribución de la variable peso, podemos observar un comportamiento que no seria explicado completamente por la relación que planteada. Posteriormente re evaluaremos (con los supuestos de los errores) esta afirmación.

Estimamos el siguiente modelo: \(Peso = \beta_0 + \beta_1 Diametro\) 

model1=lm(arboles$peso~ arboles$diametro)
stargazer(model1, type = 'html')
Dependent variable:
peso
diametro 5.103***
(0.251)
Constant -9.020***
(1.413)
Observations 90
R2 0.825
Adjusted R2 0.823
Residual Std. Error 3.435 (df = 88)
F Statistic 413.961*** (df = 1; 88)
Note: p<0.1; p<0.05; p<0.01

Ahora, se debe realizar un análisis de los errores para comprobar que se cumplen los supuestos. Existen 4 principales:

\(*\) Normalidad de los errores.

\(*\) Homocedasticidad de los errores.

\(*\) Ausencia de autocorrelación de los errores.

\(*\) Media cero de los errores. 

par(mfrow=c(2,2))
plot(model1)

X=ols_test_normality(model1)
G=matrix(c(matrix(c(round(X$shapiro$p.value,3),round(ad.test(model1$residuals)$p.value,3),'', '','' ), ncol=1), matrix(c('','',round(bptest(model1)$p.value,3),round(gqtest(model1)$p.value,3),''), ncol=1), 
matrix(c('','','','',round(t.test(model1$residuals)$p.value,3)), ncol=1)), ncol=3)
colnames(G)<- c("Normalidad","Homocedasticidad","Media cero")
rownames(G)<- c("Shapiro","Anderson","Breusch-Pagan", "Goldfeld-Quandt",'T-test')
  G %>%
   kbl(booktabs = T,) %>%
   kable_classic_2(full_width = F)
Normalidad Homocedasticidad Media cero
Shapiro 0.003
Anderson 0.025
Breusch-Pagan 0.001
Goldfeld-Quandt 0.012
T-test 1

Para validar la normalidad de los errores se puede hacer uso de dos herramientas. La primera es realizar gráficas de normalidad y, lo segundo es, apoyarnos con test estadísticos, como Shapiro wilk y Anderson-Darling. Con los test (se tiene un valor p menor a 0.05, por lo cual se rechaza la hipótesis nula de normalidad) y, con los gráficos de normalidad podemos decir que no existe evidencia a favor de la normalidad de los residuos de la regresión, lo cual, se estaría violando este supuesto.

Otro supuesto escencial del modelo es que, los resisuales son iid, es decir, independiente e indeticamente distribuidos, esto quiere decir que no debe existir autocorrelacion de los errores y estos deben ser homocedasticos.

Para validad el supuesto de homocedasticidad, el cual consiste en que la varianza de los errores sea contante, se realizará pruebas estadísticas, como el de Breusch-Pagan y el de Goldfeld-Quandt, comparándolos con los gráficos de Scale-Location y Residuals vs fitted. Con los resultados obtenidos, se puede decir que, existe evidencia a favor de heterocedasticidad presente en la regresión, dado que los valores p en las pruebas son menores al valor de significancia típico de 0.05, por tanto se rechaza la hipótesis nula de homocedasticidad. En términos gráficos se puede ver que los errores siguen una pendiente (de forma cuadrática) y no cuentan con un comportamiento de nube tipico de errores homocedasticos. Con lo anterior, se estaría violando uno de los supuestos necesarios y, por lo tanto, la estimación por MCO no sería eficiente, ya que la varianza no es óptima, y para este estudio, los estadísticos t y f no se usarían tal y como están.

Por otro lado, en cuanto a la no autocorrelación de los errores, la violación de este supuesto se da principalmente en muestras de serie de tiempo y no corte transversal, dado que no es posible tener autocorrelacion seriales en los errores si las observaciones se realizan todas en el mismo momento. Por ello, se asume que este supuesto se cumple y los errores no estarian correlacionados.

Adicionalmente, la media cero de los errores hace referencia a la especificación correcta del modelo, por ejemplo, si estamos ante un caso de variables con tendencia temporal (sin transformar), es de esperarse que el término error contenga el comportamiento de esta variable o el de la correcta especificación y genere un efecto constante en los residuales. Lo anterior se puede validar de dos maneras: con la gráfica Residuals vs Fitted, y, la esperanza matemática de los errores. Sin embargo en este problema, dado que estamos mirando un corte transversal, con variables sin dependencia temporal, podemos ver que la media de los errores parece ser cero, y, al usar un test estdistico, T-test, corroboramos efectivamente esta afirmación.

Por último, para corregir las diversas violaciones a los supuestos (normalidad, y homocedasticidad), se realiza lasas transformaciones de los datos, de tal manera que, la relación entre las variables (transformadas) puedan capturarse de una mejor manera. 

bc=boxcox(lm(arboles$peso~ arboles$diametro), lambda = -1:1)

(lambda <- bc$x[which.max(bc$y)])
## [1] 0.07070707
model2 <- lm(log(arboles$peso)~ arboles$diametro)

Usando la metodología \(Box-Cox\) se observa que, el \(\lambda\) óptimo para estos datos es cercano a \(0\) en un invervalo del \(95\%\) de confianza, lo cual, representa que, la variable dependiente en ese caso en peso, tiene una transformación logarítmica \(log(peso)=\beta_0+\beta_1Diametro\), resultando en un modelo log-lin.

Con esta nueva especificación del modelo se puede analizar el comportamiento de los errores y ver si se validan los supuestos de la regresión lineal (MCO)., de ser asi,entonces se podría realizar inferencia con los resultados obtenidos. 

par(mfrow=c(2,2))
plot(model2)

X=ols_test_normality(model2)
G=matrix(c(matrix(c(round(X$shapiro$p.value,3),round(ad.test(model2$residuals)$p.value,3),'', '','' ), ncol=1), matrix(c('','',round(bptest(model2)$p.value,3),round(gqtest(model2)$p.value,3),''), ncol=1), 
matrix(c('','','','',round(t.test(model2$residuals)$p.value,3)), ncol=1)), ncol=3)
colnames(G)<- c("Normalidad","Homocedasticidad","Media cero")
rownames(G)<- c("Shapiro","Anderson","Breusch-Pagan", "Goldfeld-Quandt",'T-test')
  G %>%
   kbl(booktabs = T,) %>%
   kable_classic_2(full_width = F)
Normalidad Homocedasticidad Media cero
Shapiro 0.334
Anderson 0.582
Breusch-Pagan 0.049
Goldfeld-Quandt 0.321
T-test 1

Con estas nuevas gráficas y análisis estadísticos se logra evidenciar una mejora significativa en los supuestos. En primera instancia, cuando se observa las gráficas de Residual vs Fitted se tiene una línea horizontal y unos puntos con comportamiento de nube, lo cual, muestra aún que los residuales son lineales y parece tener errores homocedasticos. Cuando analizamos la grafica qq de normalidad, los puntos estan muy cercanos a la linea de \(45\) grados mostrando que, los errores son normales.

Cuando se analiza el scale-location, se corrobora lo visto en la gráfica de Residual vs Fitted, donde parece no existir heterocedasticidad en los errores.

Estos comportamientos en las gráficas están también soportados con resultados de los diversos test estadísticos presentados en la tabla, los cuales evidencian que, a diferencia de la regresión sin ninguna transformación, en ésta no se tiene problemas de no normalidad, ni de heterocedasticidad y, por lo tanto, se puede decir que los coeficientes son insesgados, consistentes y eficientes y el análisis que se realizará a partir de estos resultados, el cuál tendrá las características mencionadas. 

stargazer(model1,
          model2,
          title="Estimación de peso de los árboles",
          type = "html",
          float = TRUE,
          report = "vcs*",
          no.space = TRUE,
          header=FALSE,
          single.row = FALSE,
          font.size = "small",
          intercept.bottom = F,
          column.labels = c("Normal", "Transformado"),
          column.separate = c(1, 1),
          digits = 4,
          t.auto = F,
          p.auto = F,
          notes.align = "l",
          notes.append = FALSE,
          omit.stat=c("f")
          )
Estimación de peso de los árboles
Dependent variable:
peso peso)
Normal Transformado
(1) (2)
Constant -9.0203 1.3280
(1.4129)*** (0.0598)***
diametro 5.1026 0.2782
(0.2508)*** (0.0106)***
Observations 90 90
R2 0.8247 0.8865
Adjusted R2 0.8227 0.8852
Residual Std. Error (df = 88) 3.4348 0.1453
Note: p<0.1; p<0.05; p<0.01

Como se observa, el modelo transformado tiene un mayor grado explicativo, es decir, un mayor \(R^2 Ajustado\), que el modelo sin transformar y, adicionalmente, se cumplen los supuestos de los errores.

En el modelo transformado los coeficientes de la constante y el diámetro son significativos al \(95\%\) de confianza, la constante \(\beta_0\) es de \(1.3064\) lo cual indica que, aún si el arbol no tiene diametro, su peso tomará el valor transformado \(exp(1.3064)=3.773417\).

Por su parte, el coeficiente \(\beta_1\) que acompaña a la variable diámetro nos dice que, por cada metro adicional se incrementa el peso de los árboles en un \(32.073\%\).

El modelo planteado tiene un alto grado explicativo, si se observa el \(R^2\) y el \(R^2 Ajustado\), se explica el \(0.8871\) y \(0.8858\) respectivamente, de la variable dependiente (peso).

Luego de analizar los coeficientes del modelo planteado, se realiza un análisis de pronóstico dentro de muestra. Se toma un sub set de la base de datos \(Arboles\) y se genera una validación cruzada, para evaluar el poder predictivo del modelo. 

p=predict(model2,list(diametro=8),interval = "confidence",level = 0.95)

G=matrix(c(p[1], p[2], p[3]), ncol=3)
colnames(G)<- c("Media","Inferior","Superior")
rownames(G)<- c("Peso estimado en porcentaje")

  G %>%
   kbl(booktabs = T,) %>%
   kable_classic_2(full_width = F)
Media Inferior Superior
Peso estimado en porcentaje 2.635448 2.802359 2.663267

Como se evidencia en la gráfica y en el MAPE, (mean absolut percent error), el modelo tiene algo de poder predictivo aunque no es del 100 porciento, con este modelo transformado, es posible predecir el \(88.56\%\) de los datos reales del peso de los árboles, resultando un error de pronóstico del \(11.43\%\). Por último, aunque éste no es el mejor en términos predictivos( tal vez por que es un modelo univariado) si tiene un buen ajuste y cumple los supuestos principales de los errores bajo las estimaciones de MCO, generando estimadores consistentes, insesgados y eficientes.