Entregable

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.

Metodologia

Para el desarrollo de este caso de estudio se seguirán los siguientes pasos:

Descripción de las variables estudiadas.
Definición y Calculo de modelos y gráficas para validación de supuestos.
Generación de matriz de resultados e hipótesis para validación de supuestos.
Evaluación de modelos y supuestos.
Definición e interpretación del modelo seleccionado.

Para la evaluación de los modelos se desarrolará la siguiente metodologia:

Validación de supuestos del error
Validación de linealidad del modelo
Validación de significancia de los coeficientes

Para los modelos que pasen las validaciones se determinará la escogecia del mejor modelo utilizando el R2.

1. Descripción de las variables estudiadas.

1.1 Estadística descriptiva:

library(paqueteMOD)
data(arboles)

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

Interpretación: El data set a analizar cuenta con 90 registros; con respecto a la variable dependiente (peso) se tiene una media de 18.77 con una distribución heterogénea (CV > 0.3) por otro lado, las variables independientes, altura y diámetro, presentan valores medios de 6.63 y 5.45 respectivamente, con distribuciones homogéneas, ambas con CV < 0.3.

1.2 Relacion entre la variable dependiente y posibles predictores:

library(patchwork)
library(ggplot2)

g1 = ggplot(arboles, aes(x=altura, y=peso)) + geom_point()
g2 = ggplot(arboles, aes(x=diametro, y=peso)) + geom_point()

g1 + g2

Interpretación: Con respecto a la relación de Peso vs Altura se identifica gráficamente una correlación posiblemente lineal y positiva, por otro lado, con respecto a la relación de Peso Vs Diámetro a pesar de no descartar una posible relación lineal y positiva, se considera que los datos refleja una posible relación cuadrática.

1.3 Distribucíón de las variables

g1 = ggplot(arboles, aes(x=altura)) + geom_boxplot()+ coord_flip()
g2 = ggplot(arboles, aes(x=diametro)) + geom_boxplot()+ coord_flip()
g3 = ggplot(arboles, aes(x=peso)) + geom_boxplot()+ coord_flip()

g4 =  ggplot(arboles) + 
      geom_histogram(aes(x=altura), bins = 10, alpha = 0.5, color = 'black') +
      labs(x = 'Altura', y = 'Cuenta')

g5 =  ggplot(arboles) + 
      geom_histogram(aes(x=diametro), bins = 10, alpha = 0.5, color = 'black') +
      labs(x = 'Diametro', y = 'Cuenta')

g6 =  ggplot(arboles) + 
      geom_histogram(aes(x=peso), bins = 10, alpha = 0.5, color = 'black') +
      labs(x = 'Peso', y = 'Cuenta')

(g1 + g2 + g3)/(g4 + g5 + g6)

Interpretación: Con respecto a la distribución de las variables, la Altura presenta una distribución ligeramente sesgada a la derecha sin presencia de outliers, el Diámetro presenta una distribución aproximadamente normal sin presencia de outliers y por ultimo, la variable dependiente (Peso) presenta una distribución aproximadamente normal y contiene algunos pocos datos atípicos.

1.4 Normalidad en las variables

library(pander)

pander(shapiro.test(arboles$altura))

Shapiro-Wilk normality test: `arboles$altura`
Test statistic	P value
0.9791	0.156

pander(shapiro.test(arboles$diametro))

Shapiro-Wilk normality test: `arboles$diametro`
Test statistic	P value
0.9911	0.8079

pander(shapiro.test(arboles$peso))

Shapiro-Wilk normality test: `arboles$peso`
Test statistic	P value
0.9278	9.028e-05 * * *

Interpretación: _ Se ejecuta la prueba de hipótesis Ho: La variable distribuye normal vs H1: La variable no distribuye normal, para cada una de las variables utilizando la prueba shapiro-Wilk; con un nivel de significancia de 0.05 no se rechaza la hipótesis nula de que las variables Arboles y Diámetro distribuyen normal pero si hay evidencia para rechazar la normalidad en la variable dependiente Peso, con el fin de encontrar normalidad en la vaiable dependiente se genera transformación logaritmica en esta: _

shapiro.test(log(arboles$peso))

## 
##  Shapiro-Wilk normality test
## 
## data:  log(arboles$peso)
## W = 0.99033, p-value = 0.7545

Interpretación: Con un nivel de significance de 0.05 nbno hay evidencia para rechazar la normalidad en la variable Log(peso )

pander(cor.test(arboles$peso, arboles$altura), method =c("spearman"))

Pearson’s product-moment correlation: `arboles$peso` and `arboles$altura`
Test statistic	df	P value	Alternative hypothesis	cor
15.68	88	3.213e-27 * * *	two.sided	0.8582

pander(cor.test(arboles$peso, arboles$diametro), method =c("spearman"))

Pearson’s product-moment correlation: `arboles$peso` and `arboles$diametro`
Test statistic	df	P value	Alternative hypothesis	cor
20.35	88	4.978e-35 * * *	two.sided	0.9081

2. Definición y Calculo de modelos y gráficas para validación de supuestos.

Se desarrollaran los siguientes modelos para ser evaluados y comparados en su desempeño bajo el metodo de Minimos Cuadraros:

Modelo 1: Log Peso vs Altura Modelo 2: Peso vs Altura Modelo 3: Log Peso vs Diametro Modelo 4: peso vs Diametro

2.1 Modelo Log Peso vs Altura (Modelo 1)

library(pander)

Modelo1 =  lm(log(peso) ~ altura, arboles)

pander(summary(Modelo1))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	1.491	0.09055	16.46	1.321e-28
altura	0.2038	0.01318	15.47	7.895e-27

Fitting linear model: log(peso) ~ altura
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
90	0.2237	0.7311	0.728

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

library(ggplot2)
ggplot(arboles, aes(altura, log(peso))) + 
         geom_point() + 
         geom_smooth(method = "lm", level = 0.95, formula = y ~ x)

Interpretación Inicial: De un análisis visual se podría pensar que los errores tienen media 0, no parecen tener varianza constante y distribuyen normal, estos supuestos se validaran con pruebas de hipótesis en la siguiente sección. Por otro lado se puede apoyar la linealidad del de la relación

2.2 Modelo Peso vs Altura (Modelo 2)

library(pander)

Modelo2 =  lm(peso ~ altura, arboles)

pander(summary(Modelo2))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-7.046	1.705	-4.133	8.136e-05
altura	3.891	0.2481	15.68	3.213e-27

Fitting linear model: peso ~ altura
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
90	4.211	0.7365	0.7335

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

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

Interpretación Inicial: De un análisis visual se podría pensar que los errores no tienen media 0 (Se nota una mayor acumulación de errores negativos), no parecen tener varianza constante (En la grafica de Residuals vs Leverage se evidencia un efecto megáfono), distribuyen normal, desde luego estos supuestos se validaran con pruebas de hipótesis en la siguiente sección.

2.3 Modelo Log Peso vs Diametro (Modelo 3)

library(pander)

Modelo3 =  lm(log(peso) ~ diametro, arboles)

pander(summary(Modelo3))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	1.328	0.05977	22.22	7.488e-38
diametro	0.2782	0.01061	26.22	2.334e-43

Fitting linear model: log(peso) ~ diametro
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
90	0.1453	0.8865	0.8852

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

library(ggplot2)
ggplot(arboles, aes(altura, log(peso))) + 
         geom_point() + 
         geom_smooth(method = "lm", level = 0.95, formula = y ~ x)

Interpretación Inicial: De un análisis visual se podría pensar que los errores tienen media 0, varianza constante y distribuyen normal, estos supuestos se validaran con pruebas de hipótesis en la siguiente sección.

2.4 Modelo Peso vs Diametro (Modelo 4)

library(pander)

Modelo4 =  lm(peso ~ diametro, arboles)

pander(summary(Modelo4))

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-9.02	1.413	-6.384	7.857e-09
diametro	5.103	0.2508	20.35	4.978e-35

Fitting linear model: peso ~ diametro
Observations	Residual Std. Error	$R^2$	Adjusted $R^2$
90	3.435	0.8247	0.8227

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

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

Interpretación Inicial: De un análisis visual se podría pensar que los errores tienen media 0, no parecen tener independencia entre ellos, parecen tener varianza constante y distribuyen normal, estos supuestos se validaran con pruebas de hipótesis en la siguiente sección.

3. Generación de matriz de resultados e hipótesis para validación de supuestos

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

Modelos= c(
          "Variables : ",
          "Ho: x ~ Normal",
          "Ho: y ~ Normal ",
          "Ho: e ~ Normal ",
          "H1: Media de error <> 0 ",
          "H1: Autocorrelacion e > 0 ",
          "H1: Varianza NO constante ",
          "Linealidad del Modelo ")

Modelo_1= c(
          "Log(peso) vs altura",
          round(shapiro.test(arboles$altura)$p,4),
          round(shapiro.test(log(arboles$peso))$p,4),
          round(shapiro.test(Modelo1$residuals)$p,4),
          round(t.test(Modelo1$residuals, mu = 0)$p.value,4),
          round(dwtest(Modelo1)$p,4),
          round(gqtest(Modelo1)$p.value,4),
          "OK ")

Modelo_2= c(
          "peso vs altura",
          round(shapiro.test(arboles$altura)$p,4),
          round(shapiro.test(arboles$peso)$p,4),
          round(shapiro.test(Modelo2$residuals)$p,4),
          round(t.test(Modelo2$residuals, mu = 0)$p.value,4),
          round(dwtest(Modelo2)$p,4),
          round(gqtest(Modelo2)$p.value,4),
          "OK ")

Modelo_3= c(
          "Log(peso) vs diametro ",
          round(shapiro.test(arboles$diametro)$p,4),
          round(shapiro.test(log(arboles$peso))$p,4),
          round(shapiro.test(Modelo3$residuals)$p,4),
          round(t.test(Modelo3$residuals, mu = 0)$p.value,4),
          round(dwtest(Modelo3)$p,4),
          round(gqtest(Modelo3)$p.value,4),
          "OK ")

Modelo_4= c(
          "peso vs diametro",
          round(shapiro.test(arboles$diametro)$p,4),
          round(shapiro.test(arboles$peso)$p,4),
          round(shapiro.test(Modelo4$residuals)$p,4),
          round(t.test(Modelo4$residuals, mu = 0)$p.value,4),
          round(dwtest(Modelo4)$p,4),
          round(gqtest(Modelo4)$p.value,4),
          "OK")




results =  data.frame(Modelos, Modelo_1, Modelo_2, Modelo_3, Modelo_4)

library(kableExtra)

  kbl(results) %>%
    
  kable_material(c("striped","hover"))

Modelos	Modelo_1	Modelo_2	Modelo_3	Modelo_4
Variables :	Log(peso) vs altura	peso vs altura	Log(peso) vs diametro	peso vs diametro
Ho: x ~ Normal	0.156	0.156	0.8079	0.8079
Ho: y ~ Normal	0.7545	1e-04	0.7545	1e-04
Ho: e ~ Normal	0.0168	0.0032	0.3338	0.0028
H1: Media de error <> 0	1	1	1	1
H1: Autocorrelacion e > 0	0	0	0	0
H1: Varianza NO constante	0.9972	0.3761	0.3206	0.012
Linealidad del Modelo	OK	OK	OK	OK

   | Epic Games  | M             | 0            |

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer

stargazer(Modelo1, Modelo2, Modelo3, Modelo4,  type="text", df=FALSE)

## 
## ===============================================================
##                                 Dependent variable:            
##                     -------------------------------------------
##                     log(peso)     peso    log(peso)     peso   
##                        (1)        (2)        (3)        (4)    
## ---------------------------------------------------------------
## altura               0.204***   3.891***                       
##                      (0.013)    (0.248)                        
##                                                                
## diametro                                   0.278***   5.103*** 
##                                            (0.011)    (0.251)  
##                                                                
## Constant             1.491***  -7.046***   1.328***  -9.020*** 
##                      (0.091)    (1.705)    (0.060)    (1.413)  
##                                                                
## ---------------------------------------------------------------
## Observations            90         90         90         90    
## R2                    0.731      0.737      0.887      0.825   
## Adjusted R2           0.728      0.734      0.885      0.823   
## Residual Std. Error   0.224      4.211      0.145      3.435   
## F Statistic         239.249*** 245.977*** 687.562*** 413.961***
## ===============================================================
## Note:                               *p<0.1; **p<0.05; ***p<0.01

4. Evaluación de modelos y supuestos.

Para la evaluación de los modelos se validaran los supuestos del error y significancia de los coeficientes bajos las siguientes hipótesis y tests, adicional a esto se determinará si el modelo es valido para ser seleccionado:

Ho: e ~ Normal vs H1: e ~ No Normal | Prueba usada: Shapiro test Ho: Media e = 0 Vs H1: Media e <> o | Prueba usada: t test Ho: Autocorrelación = 0 vs H1: Autocorrelación > 0 | Prueba usada: Durbin Watson Ho: Sigma constante vs H1: Sigma No constate | Prueba usada: Goldelf Quandt

Ho: Beta0 = 0 | Beta0 <> 0 | Prueba usada: Intervalo de confianza H1: Beta1 = 0 | Beta1 <> 0 | Prueba usada: Intervalo de confianza

Para todos los tests se usa un nivel de significance del 0.05

Evaluación Modelo 1: Con un nivel de significancia de 0.05, no existe evidencia para rechazar la hipótesis de normalidad en la distribución de los errores, no existe evidencia que apoye que la media de los errores es diferente de 0, Existe evidencia para apoyarla hipótesis que los errores son independientes y no existe evidencia suficiente para rechazar la hipótesis de homocedasticidad, por lo que concluimos que este modelo es VALIDO

Evaluación Modelo 2: Con un nivel de significancia de 0.05, existe evidencia para rechazar la hipótesis de normalidad en la distribución de los errores, no existe evidencia que apoye que la media de los errores es diferente de 0, Existe evidencia para apoyarla hipótesis que los errores son independientes y no existe evidencia suficiente para rechazar la hipótesis de homocedasticidad, por lo que concluimos que este modelo es INVALIDO

Evaluación Modelo 3: Con un nivel de significancia de 0.05, no existe evidencia para rechazar la hipótesis de normalidad en la distribución de los errores, no existe evidencia que apoye que la media de los errores es diferente de 0, Existe evidencia para apoyarla hipótesis que los errores son independientes y no existe evidencia suficiente para rechazar la hipótesis de homocedasticidad, por lo que concluimos que este modelo es VALIDO

Evaluación Modelo 4: Con un nivel de significancia de 0.05, no existe evidencia para rechazar la hipótesis de normalidad en la distribución de los errores, no existe evidencia que apoye que la media de los errores es diferente de 0, Existe evidencia para apoyarla hipótesis que los errores son independientes y existe evidencia suficiente para apoyar la existencia de heterocedasticisas en los errores, por lo que concluimos que este modelo es INVALIDO

5. Definición e interpretación del modelo seleccionado.

Del punto anterior definimos que solo son validos los modelos 1 (Log peso vs Altura) y 3 (Log peso vs Diámetro), la elección del modelo se hará a través del R2 el cual nos indica que porcentaje de la variabilidad de la variable de respuesta es explicado por la variable independiente. El modelo 1 tiene un R2 de 0.731 mientras que el modelo 3 tiene un R2 de 0.887, dicho esto, el modelo seleccionado para predecir el peso de los arboles es el modelo 3.

Interpretación del modelo seleccionado:

El modelo seleccionado es Log(peso) = 1.328 + 0.278*(Diámetro)

Los coeficientes de este modelo se interpretan como:

Por cada unidad incrementada en el diámetro del árbol el logaritmo del peso se incrementa en 0.278
A pesar de que no tiene sentido el peso de un árbol con diámetro 0, se podría analizar el coeficiente independiente como 1.328 el log del peso del árbol cuando el diámetro es despreciable (muy pequeño) o muy cercano a 0.

Por otro lado, los intervalos de confianza para los parámetros de estos estimadores son:

Beta0 = 1.328 +- 0.030, esto quiere decir (1.268, 1.388) Beta1 = 0.278 +- 0.011, esto quiere decir (0.267, 0.289)

Dado que el cero no esta incluido en los intervalos de confianza se puede concluir que ambos rechazan la hipótesis de nulidad.

Por otro lado, se puede decir que nuestro modelo seleccionado explica el 88.7% de la variabilidad del log del peso a través de el diámetro.