Caso Biomasa - Métodos Estadísticos para la Toma de Decisiones

Actividad 1: Caso Biomasa

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.

library(paqueteMOD)
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.1.3

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.1.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.1.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(MASS)

## Warning: package 'MASS' was built under R version 4.1.3

library(tinytex) 

## Warning: package 'tinytex' was built under R version 4.1.3

Este modelo utilizará una base de datos de “arboles”, 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.

Carga y visualización de los datos

data("arboles")
head(arboles)

## # A tibble: 6 x 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

Altura:

• Media: La altura promedio de los árboles es de 6,63.
• Mediana: la mitad de los árboles tienen altura inferior o igual a 6.45 y la otra mitad tienen una altura superior a los 6.45.
• Desviación estándar: al altura del árbol está relativamente cerca de la media, con una variación de aproximadamente 1.80.

Diámetro:

• Media: El diámetro promedio de los árboles es de 5.45.
• Mediana: la mitad de los árboles tienen un diámetro inferior o igual a 5.40 y la otra mitad tienen una altura superior a los 5.40.
• Desviación estándar: el diámetro del árbol está relativamente cerca de la media, con una variación de aproximadamente 1.45.

Peso:

• Media: El peso promedio de los árboles es de 18.77.
• Mediana: la mitad de los árboles tienen un peso inferior o igual a 17.48 y la otra mitad tienen una altura superior a los 17.48.
• Desviación estándar: el peso del árbol está relativamente cerca de la media, con una variación de aproximadamente 8.16.

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.1.3

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:MASS':
## 
##     area

p1 = ggplot(arboles, aes(x=diametro))+geom_boxplot()
p2 = ggplot(arboles, aes(x=altura))+ geom_boxplot()
p3 = ggplot(arboles, aes(x=diametro))+geom_density()
p4 = ggplot(arboles, aes(x=altura))+ geom_density()

(p1 + p2)/(p3 + p4)

Los diagramas de caja de diámetro y altura indican que la muestra exhibe valores atípicos cuando se muestra fuera del bigote.

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

Diagrama de dispersión

library(ggplot2)

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

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

g1 + g2

Existe una relación lineal positiva entre la variable peso (y) y la variable altura (x), indicando que a mayor altura del árbol, mayor peso. La variable peso (y) y la variable diámetro (x) tienen una relación lineal positiva, lo que significa que cuanto mayor sea el diámetro del árbol, mayor será el peso.

Normalidad de las variables

shapiro.test(arboles$altura)

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

De la prueba SW a un nivel de significancia del 5%, se puede suponer que la variable altura proviene de una distribución normal debido a que el valor p es mayor a 0.05.

shapiro.test(arboles$diametro)

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

De la prueba SW a un nivel de significancia del 5%, se puede suponer que la variable diámetro proviene de una distribución normal debido a que el valor p es mayor a 0.05.

shapiro.test(arboles$peso)

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

De la prueba SW a un nivel de significancia del 5%, se puede suponer que la variable peso proviene de una distribución normal debido a que el valor p es mayor a 0.05.

Correlación de las variables

cor(arboles$altura, arboles$peso)

## [1] 0.8582009

La altura tiene una correlación de 0.85, lo que indica una relación lineal positiva fuerte frente al peso del árbol.

cor(arboles$diametro, arboles$peso)

## [1] 0.908123

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.

Inferencia sobre la correlación

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.

library('corrplot') 

## Warning: package 'corrplot' was built under R version 4.1.3

## corrplot 0.92 loaded

attach(arboles)
arboles <- arboles[,-1]
head(arboles)

## # A tibble: 6 x 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

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

correlacion<-round(cor(arboles),2)
correlacion

##          peso diametro altura
## peso     1.00     0.91   0.86
## diametro 0.91     1.00   0.94
## altura   0.86     0.94   1.00

En 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.

Estimación del modelo por Mínimos Cuadrados Ordinarios

- Diametro

modelo1=lm(peso ~ diametro, arboles)
summary(modelo1)

## 
## 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

Validación de los supuestos sobre los errores

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

Supuesto de normalidad de los errores

e= modelo1$residuals

# Ho: los errores tienen distribución normal
shapiro.test(e)

## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.95356, p-value = 0.002793

Supuesto de independencia de errores

library(lmtest)
# Ho: los errores sin independientes ( no estan relacionados)
dwtest(modelo1)

## 
##  Durbin-Watson test
## 
## data:  modelo1
## DW = 1.0719, p-value = 1.035e-06
## alternative hypothesis: true autocorrelation is greater than 0

Supuesto de varianza constante

# Ho: la varianza de los erroes es constante
gqtest(modelo1)

## 
##  Goldfeld-Quandt test
## 
## data:  modelo1
## GQ = 2.0131, df1 = 43, df2 = 43, p-value = 0.01196
## alternative hypothesis: variance increases from segment 1 to 2

- Altura

modelo2=lm(peso ~ altura, arboles)
summary(modelo2)

## 
## 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

Validación de los supuestos sobre los errores

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

Supuesto de normalidad de los errores

e= modelo2$residuals

# Ho: los errores tienen distribución normal
shapiro.test(e)

## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.95439, p-value = 0.003157

Supuesto de independencia de errores

library(lmtest)
# Ho: los errores sin independientes ( no estan relacionados)
dwtest(modelo2)

## 
##  Durbin-Watson test
## 
## data:  modelo2
## DW = 0.90216, p-value = 4.234e-09
## alternative hypothesis: true autocorrelation is greater than 0

Supuesto de varianza constante

# Ho: la varianza de los erroes es constante
gqtest(modelo2)

## 
##  Goldfeld-Quandt test
## 
## data:  modelo2
## GQ = 1.1017, df1 = 43, df2 = 43, p-value = 0.3761
## alternative hypothesis: variance increases from segment 1 to 2

Supuesto de linealidad del modelo

- Representación gráfica del modelo

gr1 <- ggplot(arboles, aes(diametro, peso)) +
            geom_point() + 
            geom_smooth(method = "lm",  level = 0.95, formula = y ~ x)

gr2 <- ggplot(arboles, aes(altura, peso)) +
            geom_point() + 
            geom_smooth(method = "lm",  level = 0.95, formula = y ~ x)

gr1 + gr2

# Tranformaciones

- Modelo log-lin

modelo3=lm(log(peso) ~ diametro, arboles)
summary(modelo3)

## 
## Call:
## lm(formula = log(peso) ~ 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 ***
## 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

- Modelo lin-log

modelo4=lm(peso ~ log(diametro), arboles)
summary(modelo4)

## 
## Call:
## lm(formula = peso ~ log(diametro), data = arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6316 -3.3655 -0.6562  2.2549 16.9576 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -19.909      2.629  -7.573 3.43e-11 ***
## log(diametro)   23.369      1.564  14.941  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.362 on 88 degrees of freedom
## Multiple R-squared:  0.7172, Adjusted R-squared:  0.714 
## F-statistic: 223.2 on 1 and 88 DF,  p-value: < 2.2e-16

- Modelo log-log

modelo5=lm(log(peso) ~ log(diametro), arboles)
summary(modelo5)

## 
## Call:
## lm(formula = log(peso) ~ log(diametro), data = arboles)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3150 -0.1115 -0.0170  0.1064  0.3327 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.61846    0.09793   6.315 1.07e-08 ***
## log(diametro)  1.34405    0.05826  23.070  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1625 on 88 degrees of freedom
## Multiple R-squared:  0.8581, Adjusted R-squared:  0.8565 
## F-statistic: 532.2 on 1 and 88 DF,  p-value: < 2.2e-16

Comparación de modelos

Al comparar los valores de \(R^2\), el modelo3 presenta el mayor porcentaje de explicación de la variabilidad de Y.

Evaluación de la linealidad del modelo

library(paqueteMOD)
data(ventas) 
modelo3=lm(log(peso) ~ diametro, arboles)
summary(modelo3)

## 
## Call:
## lm(formula = log(peso) ~ 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 ***
## 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.1453, menor a la de los anteriores modelos, y el coeficiente de determinación $ R^2 = 0.8865$ indicando que el diámetro explica el 88.6% de las variaciones totales del logaritmo del peso.

par(mfrow = c(1,2))
bc1<-boxcox(lm(log(peso) ~ diametro, arboles), lambda = -2:2)
#Se repite el proceso pero esta vez estrechando el rango de valores de lambda 
bc2<-boxcox(lm(log(peso) ~ diametro, arboles), lambda = 0:2)

(lambda <- bc2$x[which.max(bc2$y)])

## [1] 1.252525

bptest(modelo3)$p.value

##         BP 
## 0.04889243

La primera gráfica muestra el comportamiento aleatorio de los residuos en relación con los valores ajustados, proporcionando evidencia gráfica para respaldar la suposición de varianza constante. Por otro lado, el p-valor de la prueba de compensación de Breusch es de 0.048, lo que indica que no se cumple el supuesto al 5% de significancia, pero su p-valor es muy cercano.

shapiro.test(modelo3$residuals)$p.value

## [1] 0.3338082

El segundo gráfico muestra que, excepto por algunos valores más bajos, los residuos se ubican mejor en la línea de valores teóricos. La prueba de Shapiro también da un valor de 0,333, que es superior a 0,05, por lo que la significancia es del 5%, lo que no rechaza el supuesto de normalidad del residual.

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

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