1 Introduccion

Los modelos de estimación de biomasa pueden influenciar estimaciones locales, regionales y globales de los volúmenes de madera de los bosques de acuerdo a su robustez y precisión. Los estándares actuales de procedimientos operacionales para la medición de carbono recomiendan el uso de modelos existentes (validados por expertos). Sin embargo, es difícil de conocer qué tan reproducibles pueden ser los modelos debido a ciertas variedades de climas y microclimas de los bosques.

En este estudio se probaran diferentes combinaciones de las variables predictoras (medidas de los arboles llamadas variables dendrometricas) haciendo uso de la regresion lineal por minimos cuadrado ordinarios para estimar el peso de los arboles (AGB, above ground biomass).

Se uso la metrica del Error Cuadratico Medio (RMSE, por sus siglas en ingles) para comparar el error entre las estimaciones de los modelos para determinar cual es el mas preciso.

2 Paquetes de R y Datos

2.1 Paquetes

library(dplyr)
library(readr)
library(corrplot)
library(rstatix)
library(ggplot2)
library(ggpubr)
library(dplyr)
library(Metrics)

dplyr: para manipular las tablas.
readr: para leer los datos en R.
corrplot: para mostrar la matriz de correlacion de una manera comprensible.
rstatix: pruebas estadisticas parametricas.
ggplot2: para graficar las distribuciones y las curvas de regresión de los datos.
ggpubr: para tests estadísticos visuales.
Metrics: para calcular el Error Cuadratico Medio de los modelos.

2.2 Datos

trees <- read_csv("../Datos/trees.csv")
trees

arbol: valor de indice de los arboles.
diam.bas: diametro basal del arbol medido a $30cm$ sobre el nivel del piso en $cm$.
altura: altura total del arbol en $metros$.
area.copa: area de cobertura de la copa del arbol en $m^3$ considerando que es un ovalo.
anilllos: edad del arbol, medida proveniente del conteo de anillos de un disco basal de cada arbol.
pesos: peso total en $kg$ de cada arbol (tambien refernciado como AGB).

2.3 Tema de ggplot2

# custom ggplot theme
require(extrafont, quietly = TRUE) # Fuente Poppins necesaria
loadfonts()
custom_theme <- theme(
  axis.text = element_text(size = 15),
  axis.title = element_text(size = 16),
  legend.text = element_text(size = 15),
  legend.title = element_text(size = 15),
  legend.position = "none",
  panel.background = element_blank(),
  panel.border = element_rect(size = 1, fill = NA),
  plot.title = element_text(size = 15, face = "bold", hjust = 0.5),
  text = element_text(family = "Poppins") # Fuente Poppins
)

3 Analisis Exploratorio de Datos

EDA (por sus siglas en ingles) es un proceso en el que calculamos y analizamos estadisticas sobre las variables, realizamos figuras para encontrar tendencias, relaciones, anomalias, patrones o relaciones entre los datos. El objetivo del EDA es intentar entender que nos quieren decir nuestros datos.

3.1 Correlacion

trees_cor <- cor(select(trees, -arbol)) # diagrama de correlacion sin indice
corrplot(trees_cor, method = "number", type = "upper", tl.col = "black", 
         tl.srt = 0)

El diagrama indica que hay una alta correlacion positiva entre diametro basal, altura, area de copa y el peso de los arboles. Y tambien una casi nula relacion entre la edad del arbol con su peso (esto puede deberse a las diferentes condiciones del terreno donde se ubicaban los arboles, como disponibilidad a agua y nutrientes).

tibble(
  Correlaciones = c("Basal diameter ~ AGB", "Height ~ AGB",
                    "Cup area ~ AGB", "Tree rings ~ AGB"),
  Pearson = c(cor_test(trees, diam.bas, pesos)$cor,
              cor_test(trees, altura, pesos)$cor,
              cor_test(trees, area.copa, pesos)$cor,
              cor_test(trees, anillos, pesos)$cor),
  Valor_p = c(cor_test(trees, diam.bas, pesos)$p,
              cor_test(trees, altura, pesos)$p,
              cor_test(trees, area.copa, pesos)$p,
              cor_test(trees, anillos, pesos)$p)
  )

Los valores p de las correlaciones entre diametro basal, altura y area de copa con el peso de los arboles son significativos ($p<\alpha$, $\alpha=0.05$). Podemos concluir que las variables estan altamente correlacionadas con el peso de los arboles, excepto la edad, la cual su valor $\rho$ es casi 0.

Mirando estos valores, es recomendable usar solo las variables diametro basal, altura y area de copa para explicar el peso de los arboles en los modelos de regresion.

3.2 Distribuciones y relaciones entre variables

3.2.1 Diametro basal

plot_diam1 <- trees %>% 
  ggplot(aes(x = diam.bas)) +
  geom_histogram(bins = 10, fill = "white", col = "black") +
  labs(x = "Diam. Bas (cm)", y = "N° de arboles") +
  custom_theme

plot_diam2 <- trees %>% 
  ggplot(aes(x = diam.bas, y = pesos)) +
  geom_point(size = 3) +
  labs(x = "Diam. Bas (cm)", y = "AGB (Kg)") +
  custom_theme

ggarrange(plot_diam1, plot_diam2)

Hay una mayor cantidad de arboles con diametro basal menor a los 40 cm, de los cuales su peso sobre el nivel del suelo es menor que 1000 kg. El diagrama de puntos indica que existe una relaciones potencial entre las 2 variables, lo que en teoria modelos de biomasa es: $AGB=aDiam^b$ (donde AGB es el peso de la biomasa sobre el nivel del suelo, por sus siglas en ingles).

3.2.2 Altura

plot_altu1 <- trees %>% 
  ggplot(aes(x = altura)) +
  geom_histogram(bins = 10, fill = "white", col = "black") +
  labs(x = "Altura (m)", y = "N° de arboles") +
  custom_theme

plot_altu2 <- trees %>% 
  ggplot(aes(x = altura, y = pesos)) +
  geom_point(size = 3) +
  labs(x = "Altura (m)", y = "AGB (Kg)") +
  custom_theme

ggarrange(plot_altu1, plot_altu2)

La distribucion de la altura de los arboles parece ser bimodal, teniendo una alta frecuencia alrededor de los 8 m y otra por los 14 m. La relacion entre la altura de los arboles y su peso tambien es potencial: $AGB=aAltura^b$.

3.2.3 Area de copa

plot_area_copa1 <- trees %>% 
  ggplot(aes(x = area.copa)) +
  geom_histogram(bins = 10, fill = "white", col = "black") +
  labs(x = expression(paste("Area de copa (", m^2, ")")), 
       y = "N° de arboles") +
  custom_theme

plot_area_copa2 <- trees %>% 
  ggplot(aes(x = area.copa, y = pesos)) +
  geom_point(size = 3) +
  labs(x = expression(paste("Area de copa (", m^2, ")")), y = "Peso (Kg)") +
  custom_theme

ggarrange(plot_area_copa1, plot_area_copa2)

Se observa una gran mayoria de arboles con un area de copa de casi $50m^2$. La relacion de area de copa con el peso de los arboles tambien es potencial: $AGB=aAreaCopa^b$.

3.3 Test de normalidad

trees %>% 
  shapiro_test(diam.bas, altura, area.copa, pesos)

Segun la prueba de normalidad de Shapiro-Wilk, ninguna de las variables son significativamente normales.

qq_diam <- ggqqplot(trees, x = "diam.bas")
qq_altu <- ggqqplot(trees, x = "altura")
qq_area <- ggqqplot(trees, x = "area.copa")
qq_peso <- ggqqplot(trees, x = "pesos")

ggarrange(qq_diam, qq_altu, qq_area, qq_peso)

Los graficos qq corroboran los resultados de las pruebas de Shapiro-Wilk de forma visual. Los datos al final se separan de la recta, denotando que es necesario una transformación para obtener la linearidad.

4 Modelamiento

El objetivo de este estudio es encontrar el mejor modelo que explique el peso de los arboles (Prosopis Sp.) con mejor precision, usando la metrica RMSE (raiz del error cuadratico medio, por sus siglas en ingles) para encontrar el modelo con menor error.

4.1 Regresion lineal por minimos cuadrados

Como se vio en 3.2. Distribuciones y relaciones entre variables, existe una relacion potencial entre las variables dendrometricas y AGB. Esta relacion puede ser explicada por la ecuacion $y=b_0x^{b_1}$, y haciendo una transformacion logaritmica logramos la expresion de forma lineal: \[log(y)=log(b_0)+b_1log(x)\].

De lo cual R va a estimar el modelo univariado: \[y=b_0+b_1x\]

Y el modelo multivariado seria: \[y=b_0+b_1x_1+b_2x_2+…+b_kx_k\]

Siendo k el numero de variables predictivas.

Cabe mencionar que el $y$ estimado es en realidad $log(y)$, entonces haciendo: $e^{log(y)}$ obtendremos los verdaderos valores predichos en la escala correcta y no los valores trasformados con el logaritmo.

4.1.1 Reg. Univariante

Diametro Basal

Forma potencial: \[AGB=b_0(bas.diam)^{b_1}\]

Forma lineal: \[log(AGB)=log(b_0)+b_1log(bas.diam)\]

model1 <- lm(log(pesos) ~ log(diam.bas), data = trees)
summary(model1)

## 
## Call:
## lm(formula = log(pesos) ~ log(diam.bas), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9408 -0.1697  0.0435  0.2297  1.2847 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -2.5354     0.5776  -4.389 0.000147 ***
## log(diam.bas)   2.4979     0.1759  14.199 2.55e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4278 on 28 degrees of freedom
## Multiple R-squared:  0.8781, Adjusted R-squared:  0.8737 
## F-statistic: 201.6 on 1 and 28 DF,  p-value: 2.554e-14

Los coeficientes estimados del modelo son el Intercepto = -2.5354 y la pendiente = 2.4979, pero que en realidad son $log(b_0)$ y $b_1$ respectivamente. Con lo que sera necesaria una transformacion al primer valor $e^{log(b_0)}=b_0$ para obtener su estimado del coeficiente en la escala correcta.

data.frame(coef = c(exp(coef(model1))[1],
                    coef(model1)[2]),
           row.names = c("b0","b1"))

Con los coeficientes estimados, se puede plantear la ecuacion potencial que explica el peso de los arboles (AGB) usando su diametro basal de esta forma:

\[AGB=0.079(bas.diam)^{2.498}\] ##### RMSE {-}

# predictions
pred_model1 <- exp(predict(model1))

# RMSE
rmse_model1 <- rmse(trees$pesos, pred_model1)
data.frame(RMSE = rmse_model1, row.names = "AGB ~ bas.diam")

El RMSE (Error Cuadratico Medio) del primer modelo es 330.076

Grafico de ajuste

trees |> 
  mutate(pred_model1) |> 
  ggplot(aes(x = diam.bas)) +
  geom_point(aes(y = pesos), size = 3) +
  geom_line(aes(y = pred_model1), color = "red", size = 1.5, alpha = 0.75) +
  geom_text(label = "RMSE = 330.0760", x = 25, y = 3600) +
  labs(title = "AGB ~ Diam. Bas", 
       x = "Diam.Bas (cm)",
       y = "AGB (Kg)") +
  custom_theme

Altura

Forma potencial: \[AGB=b_0(altura)^{b_1}\]

Forma lineal: \[log(AGB)=log(b_0)+b_1log(altura)\]

model2 <- lm(log(pesos) ~ log(altura), data = trees)
summary(model2)

## 
## Call:
## lm(formula = log(pesos) ~ log(altura), data = trees)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09020 -0.34016  0.04133  0.41138  1.07915 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.3934     0.8744  -3.881 0.000578 ***
## log(altura)   4.0380     0.3903  10.345 4.53e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5579 on 28 degrees of freedom
## Multiple R-squared:  0.7926, Adjusted R-squared:  0.7852 
## F-statistic:   107 on 1 and 28 DF,  p-value: 4.526e-11

Los coeficientes estimados del modelo son el Intercepto = -3.3934 y la pendiente = 4.0380, pero que en realidad son $log(b_0)$ y $b_1$ respectivamente. Con lo que sera necesaria una transformacion al primer valor $e^{log(b_0)}=b_0$ para obtener su estimado del coeficiente en la escala correcta.

data.frame(coef = c(exp(coef(model2))[1],
                    coef(model2)[2]),
           row.names = c("b0","b1"))

Con los coeficientes estimados, se puede plantear la ecuacion potencial que explica el peso de los arboles (AGB) usando su diametro basal de esta forma:

\[AGB=0.034(bas.diam)^{4.038}\] ##### RMSE {-}

# predictions
pred_model2 <- exp(predict(model2))

# RMSE
rmse_model2 <- rmse(trees$pesos, pred_model2)
data.frame(RMSE = rmse_model2, row.names = "AGB ~ altura")

El RMSE (Error Cuadratico Medio) del segundo modelo es 584.98

Grafico de ajuste

trees |> 
  mutate(pred_model2) |> 
  ggplot(aes(x = altura)) +
  geom_point(aes(y = pesos), size = 3) +
  geom_line(aes(y = pred_model2), color = "red", size = 1.5, alpha = 0.75) +
  geom_text(label = "RMSE = 584.9777", x = 8.5, y = 3600) +
  labs(title = "AGB ~ Altura", 
       x = "Altura (m)",
       y = "AGB (Kg)") +
  custom_theme

Area de copa

Forma potencial: \[AGB=b_0(area.copa)^{b_1}\]

Forma lineal: \[log(AGB)=log(b_0)+b_1log(area.copa)\]

model3 <- lm(log(pesos) ~ log(area.copa), data = trees)
summary(model3)

## 
## Call:
## lm(formula = log(pesos) ~ log(area.copa), data = trees)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.32990 -0.49456  0.02218  0.56175  1.21648 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.2252     0.7045   0.320    0.752    
## log(area.copa)   1.4388     0.1859   7.742 1.96e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6913 on 28 degrees of freedom
## Multiple R-squared:  0.6816, Adjusted R-squared:  0.6702 
## F-statistic: 59.94 on 1 and 28 DF,  p-value: 1.963e-08

Los coeficientes estimados del modelo son el Intercepto = 0.2252 y la pendiente = 1.4388, pero que en realidad son $log(b_0)$ y $b_1$ respectivamente. Con lo que sera necesaria una transformacion al primer valor $e^{log(b_0)}=b_0$ para obtener su estimado del coeficiente en la escala correcta.

data.frame(coef = c(exp(coef(model3))[1],
                    coef(model3)[2]),
           row.names = c("b0","b1"))

Con los coeficientes estimados, se puede plantear la ecuacion potencial que explica el peso de los arboles (AGB) usando su diametro basal de esta forma:

\[AGB=1.253(diam.bas)^{1.439}\] ##### RMSE {-}

# predictions
pred_model3 <- exp(predict(model3))

# RMSE
rmse_model3 <- rmse(trees$pesos, pred_model3)
data.frame(RMSE = rmse_model3, row.names = "AGB ~ area.copa")

El RMSE (Error Cuadratico Medio) del tercer modelo es 646.0653

Grafico de ajuste

trees |> 
  mutate(pred_model3) |> 
  ggplot(aes(x = area.copa)) +
  geom_point(aes(y = pesos), size = 3) +
  geom_line(aes(y = pred_model3), color = "red", size = 1.5, alpha = 0.75) +
  geom_text(label = "RMSE = 646.0653", x = 40, y = 3600) +
  labs(title = "AGB ~ Altura", 
       x = "Altura (m)",
       y = "AGB (Kg)") +
  custom_theme

4.1.2 Reg. Multivariante

De todas las combinaciones de modelos multivariantes para estos datos se procederan a mostrar los 3 mejores:

4.1.2.1 Modelo 4

Linear form: \[log(AGB)=log(b_0)+b_1log(diam.bas)+b_2log(altura)\]

model4 <- lm(log(pesos) ~ log(diam.bas) + log(altura), data = trees)
summary(model4)

## 
## Call:
## lm(formula = log(pesos) ~ log(diam.bas) + log(altura), data = trees)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78583 -0.12109  0.08681  0.19013  1.03328 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -3.1637     0.6425  -4.924 3.73e-05 ***
## log(diam.bas)   1.8639     0.3718   5.013 2.94e-05 ***
## log(altura)     1.2094     0.6326   1.912   0.0666 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4089 on 27 degrees of freedom
## Multiple R-squared:  0.8926, Adjusted R-squared:  0.8846 
## F-statistic: 112.2 on 2 and 27 DF,  p-value: 8.298e-14

Modelo estimado: \[log(AGB)=e^{-3.1637}+1.8639(diam.bas)+1.2094(altura)\] \[log(AGB)=0.0423+1.8639(diam.bas)+1.2094(altura)\]

pred_model4 <- exp(predict(model4))
rmse_model4 <- rmse(trees$pesos, pred_model4)
data.frame(RMSE = round(rmse_model4, 2), row.names = "AGB ~ diam.bas + altura")

El Error Cuadratico medio del cuarto modelo es 367.01.

4.1.2.2 Modelo 5

Forma lineal: \[log(AGB)=log(b_0)+b_1log(diam.bas)+b_2log(area.copa)\]

model5 <- lm(log(pesos) ~ log(diam.bas) + log(area.copa), data = trees)
summary(model5)

## 
## Call:
## lm(formula = log(pesos) ~ log(diam.bas) + log(area.copa), data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9794 -0.1683  0.0203  0.1938  1.2979 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.4448     0.5916  -4.133 0.000311 ***
## log(diam.bas)    2.2638     0.3366   6.726  3.2e-07 ***
## log(area.copa)   0.1799     0.2200   0.818 0.420744    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4304 on 27 degrees of freedom
## Multiple R-squared:  0.881,  Adjusted R-squared:  0.8722 
## F-statistic: 99.94 on 2 and 27 DF,  p-value: 3.311e-13

Modelo estimado: \[log(AGB)=e^{-2.4448}+2.264(diam.bas)+0.179(area.copa)\] \[log(AGB)=0.087+2.264(diam.bas)+0.179(altura)\]

pred_model5 <- exp(predict(model5))
rmse_model5 <- rmse(trees$pesos, pred_model5)
data.frame(RMSE = round(rmse_model5, 1), row.names = "AGB ~ diam.bas + area.copa")

El Error Cuadratico medio del quinto modelo es 367. Hasta ahora es el mejor modelo, ganandole al primero con RMSE = 330.1. Sin embargo, segun el ANOVA, log(area.copa) no es significativo. Dejando asi como mejor hasta ahora al primer modelo.

4.1.2.3 Modelo 6

Forma lineal: \[log(AGB)=log(b_0)+b_1log(diam.bas)+b_2log(altura)+b_3log(area.copa)\]

model6 <- lm(log(pesos) ~ log(diam.bas) + log(altura) + log(area.copa), data = trees)
summary(model6)

## 
## Call:
## lm(formula = log(pesos) ~ log(diam.bas) + log(altura) + log(area.copa), 
##     data = trees)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81334 -0.15459  0.06272  0.17729  1.00174 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -3.1418     0.6291  -4.994 3.42e-05 ***
## log(diam.bas)    1.3182     0.5191   2.540   0.0174 *  
## log(altura)      1.4713     0.6443   2.284   0.0308 *  
## log(area.copa)   0.3139     0.2129   1.474   0.1524    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4003 on 26 degrees of freedom
## Multiple R-squared:  0.9009, Adjusted R-squared:  0.8894 
## F-statistic: 78.77 on 3 and 26 DF,  p-value: 3.554e-13

Modelo estimado: \[log(AGB)=e^{-3.1418}+1.318(diam.bas)+1.471(altura)+0.314(area.copa)\] \[log(AGB)=0.043+1.318(diam.bas)+1.471(altura)+0.314(area.copa)\]

pred_model6 <- exp(predict(model6))
rmse_model6 <- rmse(trees$pesos, pred_model6)
data.frame(RMSE = round(rmse_model6, 1), row.names = "AGB ~ diam.bas + altura + area.copa")

El Error Cuadratico medio del sexto modelo es 369.7. Cabe mencionar que en el ANOVA log(area.copa) no es significativo, lo que invalidaria el modelo.

4.1.3 RMSE Modelos

# rmse data frame
data.frame(RMSE = c(rmse_model1, rmse_model2, rmse_model3,
                    rmse_model4, rmse_model5, rmse_model6),
                       row.names = c("1. AGB ~ diam.bas",
                                     "2. AGB ~ altura",
                                     "3. AGB ~ area.copa",
                                     "4. AGB ~ diam.bas + altura",
                                     "5. AGB ~ diam.bas + area.copa",
                                     "6. AGB ~ diam.bas + altura + area.copa"))

5 Conclusiones

La relacion entre las variables dendrometricas de los arboles (Prosopis Sp.) con el peso (AGB) tiene una apariencia potencial explicada por la formula: $y=b_0+b_1x$.
Es necesario realizar una transformacion logaritmica a la ecuacion: $y=b_0x^{b_1}$, con lo que obtenemos la forma lineal: $log(y)=log(b_0)+b_1log(x)$ y asi poder estimar los coeficientes de los modelos por el metodo de minimos cuadrados ordinarios.
Como hablamos de un modelo de estimacion lineal pero transformado, no es correcto usar la metrica $R^2$ o $R^2_{ajustado}$. Sin embargo, es prudente usar otra metrica como la Raiz del Error Cuadratico Medio.
El mejor modelo es el que usa la ecuacion 1.$AGB=0.079(bas.diam)^{2.498}$ que tiene menor RMSE en comparacion al resto de modelos calculados.

Estimacion de Biomasa de Prosopis Sp.

Walter Chanava

Octubre 2021

1 Introduccion

2 Paquetes de R y Datos

2.1 Paquetes

2.2 Datos

2.3 Tema de ggplot2

3 Analisis Exploratorio de Datos

3.1 Correlacion

3.2 Distribuciones y relaciones entre variables

3.2.1 Diametro basal

3.2.2 Altura

3.2.3 Area de copa

3.3 Test de normalidad

4 Modelamiento

4.1 Regresion lineal por minimos cuadrados

4.1.1 Reg. Univariante

Diametro Basal

Grafico de ajuste

Altura

Grafico de ajuste

Area de copa

Grafico de ajuste

4.1.2 Reg. Multivariante

4.1.2.1 Modelo 4

4.1.2.2 Modelo 5

4.1.2.3 Modelo 6

4.1.3 RMSE Modelos

5 Conclusiones