Caso biomasa

Liz Gutiérrez Rendón - estudiante de maestría en Ciencia de Datos

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.

#install.packages("learnr")  
#install.packages("devtools")
#install.packages("tinytex")

#devtools::install_github("dgonxalex80/paqueteMOD")
library(paqueteMOD)
library(ggplot2)
library(lmtest) # breusch pagan test 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MASS) # BOX COX
library(tinytex) 
#tinytex::install_tinytex()

1. Visualización de datos

En el presente modelo se utilizará la base de datos “árboles”, en la siguiente tabla se puede observar las variables con las que se trabajará: peso que será la variable dependiente, diámetro y altura que se serán las variables independientes. Aunque en el texto no se indica la unidad de medida de estas, vamos a suponer que el peso está medido en toneladas, el diámetro y la altura en metros.

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

2. Estadísticas descriptivas

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

Media de la altura: La altura promedio de los árboles es de 6.63 metros. Media del diámetro: El diámetro promedio de los árboles es de 5.45 metros. Media del peso: El peso promedio de los árboles es de 18.77 toneladas.

Mediana de la altura: la mitad de los árboles tienen altura inferior o igual a 6.45 metros y la otra mitad tienen una altura superior a los 6.45 metros. Mediana del diámetro: la mitad de los árboles tienen un diámetro inferior o igual a 5.40 metros y la otra mitad tienen una altura superior a los 5.40 metros. Mediana del peso: la mitad de los árboles tienen un peso inferior o igual a 17.48 toneladas y la otra mitad tienen una altura superior a los 17.48 toneladas.

Desviación estándar de la altura: al altura del árbol está relativamente cerca de la media, con una variación de aproximadamente 1.80 metros. Desviación estándar del diámetro: el diámetro del árbol está relativamente cerca de la media, con una variación de aproximadamente 1.45 metros. Desviación estándar del peso: el peso del árbol está relativamente cerca de la media, con una variación de aproximadamente 8.16 metros.

library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
p1 = ggplot(arboles, aes(x=diametro, y=peso))+geom_boxplot()
p2 = ggplot(arboles, aes(x=altura, y=peso))+ geom_boxplot()
p3 = ggplot(arboles, aes(x=diametro))+geom_density()
p4 = ggplot(arboles, aes(x=altura))+ geom_density()
(p1 + p2) / (p3 + p4)

El diagrama de cajas del diámetro y de la altura indican que la muestra tiene valores atípicos ya que estos se visualizan fuera de los bigotes.

La forma de la curva del gráfico de densidad indica que la densidad muestral fue bastante simétrica.

Normalidad de las variables

shapiro.test(arboles$diametro)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$diametro
## W = 0.99109, p-value = 0.8079

Con base en el test SW y con una significacia del 5%, se puede asumir que la variable diámetro proviene de una distribución normal porque el pvalor es mayor a 0.05.

shapiro.test(arboles$altura)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$altura
## W = 0.97909, p-value = 0.156

Con base en el test SW y con una significacia del 5%, se puede asumir que la variable altura proviene de una distribución normal porque el pvalor es mayor a 0.05.

shapiro.test(arboles$peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  arboles$peso
## W = 0.92778, p-value = 9.028e-05

Con base en el test SW y con una significacia del 5%, se puede asumir que la variable peso proviene de una distribución normal porque el pvalor es mayor a 0.05.

3. Análisis de Correlación

library(ggplot2)
gr1 <- ggplot(arboles, aes(x=altura, y=peso))+
         geom_point() + ggtitle("Altura vs Peso")

gr2 <- ggplot(arboles, aes(x=diametro, y=peso))+
         geom_point()+ ggtitle("Diametro vs Peso")

gr1 + gr2

la variable peso (y) y la variable altura (x) tienen una relación lineal positiva, lo que indica que a mayor altura del árbol mayor será su peso. La variable peso (y) y la variable diámetro (x) tienen una relación lineal positiva, lo que indica que a mayor diámetro del árbol mayor será su peso.

Correlación de las variables

cor(arboles $altura, arboles$peso)
## [1] 0.8582009
cor(arboles $diametro, arboles$peso)
## [1] 0.908123

La altura tiene una correlación de 0.85, lo que indica una relación lineal positiva fuerte frente al peso del árbol. El diámetro tiene una correlación de 0.90, lo que indica una relación lineal positiva fuerte frente al peso del árbol. La variable que mayor relación tiene con la varibale peso es el diámetro.

cor.test(arboles$altura, arboles$peso)
## 
##  Pearson's product-moment correlation
## 
## data:  arboles$altura and arboles$peso
## t = 15.684, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7918402 0.9045332
## sample estimates:
##       cor 
## 0.8582009
cor.test(arboles$diametro, arboles$peso)
## 
##  Pearson's product-moment correlation
## 
## data:  arboles$diametro and arboles$peso
## t = 20.346, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8634081 0.9386817
## sample estimates:
##      cor 
## 0.908123

Ho:ρ=0 Ha:ρ≠0 Se rechaza Ho, se concluye que ρ≠0, existe una correlación entre las variables altura y peso. Se rechaza Ho, se concluye que ρ≠0, existe una correlación entre las variables diámetro y peso.

Inferencia sobre la correlación

correlación de Pearson: El coeficiente de correlación de Pearson mide la magnitud de la asociación lineal entre dos variables numéricas en escala de razón. El resultado de esta prueba de hipótesis indica que las correlaciones entre las variables: peso, diámetro y altura son diferntes de cero.

library('corrplot') 
## corrplot 0.92 loaded
attach(arboles)
arboles <- arboles[,-1]
head(arboles)
## # A tibble: 6 × 3
##    peso diametro altura
##   <dbl>    <dbl>  <dbl>
## 1 13.7       4.7    5  
## 2 14.6       5.3    5.6
## 3 15.9       4.8    5.8
## 4  8.99      3.2    4.3
## 5  6.99      2.2    3.3
## 6 19.3       6.3    7.9

Matriz de varianzas covarianzas

var(arboles)
##              peso  diametro    altura
## peso     66.54169 10.754584 12.596798
## diametro 10.75458  2.107677  2.443919
## altura   12.59680  2.443919  3.237789

Matriz de correlación

mcor<-round(cor(arboles),2)
mcor
##          peso diametro altura
## peso     1.00     0.91   0.86
## diametro 0.91     1.00   0.94
## altura   0.86     0.94   1.00

vemos que la diagonal está conformada por unos.

install.packages("corrplot") # Instalar el paquete
## Warning: package 'corrplot' is in use and will not be installed
library(corrplot)            # Importar el paquete

corrplot.mixed(mcor)

Este gráfico muestra la matriz con los coeficientes de correlación. La diagonal son las variables, encima de ellas se encuentran circulos de colores, entre más intensidad de color mayor es la correlación, el tamaño de los círculos está asociado al valor absoluto de correlación. Por debajo de la diagonal se observan los valores exactos de correlación en colores. Lo cual para nuestro caso de estudio nos indica que la variable que mayor correlación tiene frente al peso es el diámetro con una correlación de 0.91.

Modelo de regresión lineal

gr3 <- ggplot(arboles, aes(x=diametro, y=peso)) +
     geom_point() +
     geom_smooth(method=lm) +
     theme_minimal()

gr4 <- ggplot(arboles, aes(x=altura, y=peso)) +
     geom_point() +
     geom_smooth(method=lm) +
     theme_minimal()

gr3 + gr4
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

4. Ajuste del Modelo

Estimación por mínimos cuadrados ordinarios (MCO)

## 
## Call:
## lm(formula = peso ~ diametro, data = arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3775 -2.6594  0.0237  1.8758 11.9876 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.0203     1.4129  -6.384 7.86e-09 ***
## diametro      5.1026     0.2508  20.346  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.435 on 88 degrees of freedom
## Multiple R-squared:  0.8247, Adjusted R-squared:  0.8227 
## F-statistic:   414 on 1 and 88 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = peso ~ altura, data = arboles)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.228 -1.969  0.572  2.377 15.106 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -7.0456     1.7046  -4.133 8.14e-05 ***
## altura        3.8906     0.2481  15.684  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.211 on 88 degrees of freedom
## Multiple R-squared:  0.7365, Adjusted R-squared:  0.7335 
## F-statistic:   246 on 1 and 88 DF,  p-value: < 2.2e-16

5. Comparación del modelo

Siginificancia del modelo

Con una significancia del 5%, se puede decir que los modelos son significativos para explicar a la variable peso, ya que los valores-p asociados a los coeficientes de regresión, son menores al 0.05.

El intercepto del modelo 1 es -9.02 y el del modelo 2 es -7.02. En este caso, no tienen interpretación porque las variables diámetro no toma valores de 0 (no hay árboles sin diámetro). La pendiente del modelo 1 es de 5.10 lo que significa que por cada metro de diámetro que aumente un árbol, el peso esperado aumenta en 5,10 toneladas. Para el modelo 2 el peso esperado aumentará en 3,89 toneladas.

Análisis de las métricas

El error estándar de la estimación (SEE) es una métrica importante para evaluar la calidad de ajuste del modelo. El modelo 1 arrojo un error estándar residual de 3.435, siendo menor que el error estándar del modelo 2 el cual es de 4.22.

El coeficiente de determinación \((R^2)\): es una medida de la proporción de la variabilidad de los datos explicada por la línea de regresión. En este caso, el modelo 1 arroja un \(R^2 = 0.82\), lo que indica que las variaciones del diámetro explican en un 82% de la variación total del peso de los árboles. El modelo 2 por su parte arrojó un \(R^2 = 0.736\) indicando que las variaciones de la altura explican en un 73.6% de la variabilidad total del peso de los árboles.

Hasta ahora el mejor modelo es el modelo 1 cuya variable independiente es el diámetro, puesto que arrojó valores significativos y sus medidas de desempeño fueron mayores.

6. Validación de supuestos

Verificamos gráficamente el cumplimiento de los supuestos sobre los errores.

par(mfrow = c(1,2))
plot(modelo1, which=1:2)

par(mfrow = c(1,2))
plot(modelo2, which=1:2)

Varianza Constante de los errores (Homocedasticidad)

\[ H_0: \,\, Los \,\, residuos \,\, tienen \,\, varianza \,\, constante \] \[ \,\, VS \,\, \]

\[ H_1: \,\, Los \,\, residuos \,\, NO \,\, tienen \,\, varianza \,\, constante \] Los residuales del modelo 1 parecen no cumplir con el supuesto de varianza constante, ya que tienen en el gráfico se observa un comportamiento de parábola.

Los residuales del modelo 2 parece alejarse de los valores ajustados a mayores valores, aparentando una varianza creciente. Complementariamente, con una significancia del 5% se rechaza la hipótesis de residuos con varianza constante para ambos modelos.

bptest(modelo1)$p.value
##           BP 
## 0.0006051803
bptest(modelo2)$p.value
##           BP 
## 0.0003625303

Al observar el primer gráfico

Normalidad de los errores

\[ H_0: \,\, Los \,\, residuos \,\, provienen \,\, de \,\, una \,\ distribución \,\ Normal \] \[ \,\, VS \,\, \]

\[ H_1: \,\, Los \,\, residuos \,\, NO \,\, provienen \,\, de \,\, una \,\ distribución \,\ Normal \]

Se observa que los residuales del modelo 1 tienen un ajuste no adecuado a la recta QQ plot ya que los residuales menores y mayores están por encima de la recta. En el modelo 2 los residuales menores se ubican por debaje de la recta, así como algunos residuales por encima.

shapiro.test(modelo1$residuals)$p.value
## [1] 0.002792769
shapiro.test(modelo2$residuals)$p.value
## [1] 0.003156811

En este caso, el test shapiro wilk arroja valores p menores a 0.05 por lo tanto, con una significancia del 5%, se concluye que los residuos no provienen de una distribución normal.

Independencia de los errores

En este caso no conocemos el orden en el que se recopilaron los datos, por lo tanto podríamos asumir que se hicieron bajo un marco muestral que se ejecutó respetando la independencia entre cada observación.

7. Transformación del modelo por Boxcox

Se realiza el gráfico de la curva de box cox para determinar su máximo \(\lambda_{max}\). El ajuste se hará para el modelo 1 que fue aquel con mejores métricas de exactitud.

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

(lambda <- bc$x[which.max(bc$y)])
## [1] 0.07070707
modelo3 <- lm(log(arboles$peso) ~ arboles$diametro, data = arboles)
summary(modelo3)
## 
## Call:
## lm(formula = log(arboles$peso) ~ arboles$diametro, data = arboles)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27395 -0.10180 -0.00328  0.10073  0.33742 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.32798    0.05977   22.22   <2e-16 ***
## arboles$diametro  0.27818    0.01061   26.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1453 on 88 degrees of freedom
## Multiple R-squared:  0.8865, Adjusted R-squared:  0.8852 
## F-statistic: 687.6 on 1 and 88 DF,  p-value: < 2.2e-16

La tabla de resumen, indica que el modelo 3, es significativo con una significancia del 5%. Además, se observa un error estándar es de 0.145, menor a la de los anteriores modelos, y el coeficiente de determinación $ R^2 = 0.886$ indicando que el diámetro explica el 88.6% de las variaciones totales del logaritmo del peso.

Supuestos del modelo 3 (Log(peso))

par(mfrow = c(1,2))
plot(modelo3, which = 1:2)

bptest(modelo3)$p.value
##         BP 
## 0.04889243
shapiro.test(modelo3$residuals)$p.value
## [1] 0.3338082

La primera gráfica evidencia que los residuales se comportan de manera aleatoria con respecto a los valores ajustados, dando una evidencia gráfica a favor del supuesto de varianza constante. AHora bien, el p-valor de la prueba Breusch pagan arroja un valor p de 0.048, lo que indica que no se cumple el supuesto bajo una significancia del 5%, pero su p-valor es bastante cercano.

La segunda gráfica evidencia que los residuales se localizan de una mejor manera sobre la recta de los valores teóricos, exceptuando por algunos valores menores. Además, el test de shapiro arroja un valor de 0,333 mayor que 0,05, por lo tanto con una significancia del 5%, no se rechaza la hipótesis a favor de la normalidad de los residuos.

gr4 <- ggplot(arboles, aes(x=diametro, y=log(peso))) +
     geom_point() +
     geom_smooth(method=lm) +
     theme_minimal()
print(gr4)
## `geom_smooth()` using formula = 'y ~ x'