El siguiente informe corresponde a una regresión para los datos del peso de una muestra de arboles, en función de la altura, el diametro entre otras. Los datos se muestran a continuación:

library(readxl)
datos_biomasa = read_excel("C:/Users/pacho/OneDrive/Escritorio/datos biomasa.xls")
datos_biomasa
## # A tibble: 90 x 7
##    finca   mg         peso_aerea peso_sub peso_total diametro altura
##    <chr>   <chr>           <dbl>    <dbl>      <dbl>    <dbl>  <dbl>
##  1 FINCA_1 GENOTIPO_1      12.8      0.93      13.7       4.7    5  
##  2 FINCA_1 GENOTIPO_1      13.9      0.69      14.6       5.3    5.6
##  3 FINCA_1 GENOTIPO_1      15.1      0.78      15.9       4.8    5.8
##  4 FINCA_1 GENOTIPO_1       8.08     0.91       8.99      3.2    4.3
##  5 FINCA_1 GENOTIPO_1       5.58     1.41       6.99      2.2    3.3
##  6 FINCA_1 GENOTIPO_2      18.5      0.84      19.3       6.3    7.9
##  7 FINCA_1 GENOTIPO_2      20.6      0.82      21.4       6.6    8.3
##  8 FINCA_1 GENOTIPO_2      12.7      1.08      13.8       5.3    7.3
##  9 FINCA_1 GENOTIPO_2      10.6      1.31      11.9       4.9    6.7
## 10 FINCA_1 GENOTIPO_2      15.7      0.92      16.6       5.9    7.1
## # ... with 80 more rows

Paso 1 - Plantear las variables del Modelo

Se desea evaluar la relación entre las variables peso_total del arbol como respuesta y la altura como variable predictora. Se espera en general que a mayor altura del arbol el peso de la madera del mismo se incremente.

Paso 2 - Explorar la relación entre las variables

y=datos_biomasa$peso_total
x=datos_biomasa$altura

plot(x,y)

cor(x,y)
## [1] 0.8582009

Se observa en la figura que existe una relación lineal positiva entre el peso y la altura, adicionalmente esta relación es fuerte porque el coeficiente de correlación de Pearson es de 0.85, indicando que la altura del árbol puede ser un buen predictor de su peso.

Paso 3 - Estimar el Modelo de Regresión Lineal Simple

mod_simple=lm(y~x)
mod_simple
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      -7.046        3.891

Como se observa el coeficiente beta 1 nos indica que por cada metro adicional de altura en el arbol se espera un incremento de 3.891 Tn de peso.

summary(mod_simple)
## 
## Call:
## lm(formula = y ~ x)
## 
## 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 ***
## x             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

Se observa que el modelo presenta un ajuste del 73% con base en el R2. Es decir, este modelo logra explicar el 73% de la variabilidad del peso del árbol.

Con base en el valor p se observa que la variable altura es significativa en el modelo, es decir, la altura efectivamente es una variable importante para explicar su peso.

Paso 4 - Predecir con base en el modelo

¿Se desea conocer cual es el valor estimado de la ganancia de un lote de 1000 arboles con una altura promedio de 7.2 m, si se estima que por cada Ton de madera la compañía logra producir una cantidad de productos(papel, cartón, …), cuyo valor estimado está en 60 millones de pesos (ganancia directa)?

##Escenario Promedio
y_mod=predict(mod_simple, list(x=7.2))
ganancia_arbol=y_mod*60
ganancia_total_media=ganancia_arbol*1000

#Escenario Bajo y Alto
MAE=mean(abs(mod_simple$residuals))
y_mod_min=y_mod-MAE
y_mod_max=y_mod+MAE

ganancia_total_inf=y_mod_min*60*1000
ganancia_total_sup=y_mod_max*60*1000

c(ganancia_total_inf, ganancia_total_media, ganancia_total_sup)
##       1       1       1 
## 1070983 1257986 1444990

De acuerdo con los resultados del modelo se espera que ese lote de 1000 arboles generen una ganancia de 1258 millones en un escenario medio. En un escenario bajo se esperan 1071 millones y en uno alto 1445 millones.

Modelo de Regresión Lineal Multiple

Con el objetivo de mejorar el ajuste del modelo para explicar o predecir el peso de los arboles, se incorporan otras variables predictoras adicionales, por ejemplo el diametro del árbol.

y=datos_biomasa$peso_total
x1=datos_biomasa$diametro
x2=datos_biomasa$altura

plot(datos_biomasa[,5:7])

cor(datos_biomasa[,5:7])
##            peso_total diametro    altura
## peso_total  1.0000000 0.908123 0.8582009
## diametro    0.9081230 1.000000 0.9355360
## altura      0.8582009 0.935536 1.0000000
mod_multiple1=lm(y~x1+x2)
summary(mod_multiple1)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3083 -2.5121  0.1608  2.0088 11.7446 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.1205     1.4305  -6.376 8.44e-09 ***
## x1            4.7395     0.7128   6.649 2.49e-09 ***
## x2            0.3132     0.5751   0.544    0.587    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.449 on 87 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8213 
## F-statistic: 205.5 on 2 and 87 DF,  p-value: < 2.2e-16

Como se puede observar las variables predictoras diametro y altura se encuentran altamente correlacionadas lo cual genera un problema de multicolinealidad. Esto quiere decir que en el caso del modelo multiple se observa que la variable altura no es significativa estando la variable diametro en el modelo. Esto se debe a que ya e diametro explica en parte lo que la altura ofrece al modelo (peso).

Selección de Variables

mod_multiple2=step(mod_multiple1)
## Start:  AIC=225.79
## y ~ x1 + x2
## 
##        Df Sum of Sq    RSS    AIC
## - x2    1      3.53 1038.2 224.09
## <none>              1034.7 225.79
## - x1    1    525.74 1560.5 260.76
## 
## Step:  AIC=224.09
## y ~ x1
## 
##        Df Sum of Sq    RSS    AIC
## <none>              1038.2 224.09
## - x1    1      4884 5922.2 378.80
summary(mod_multiple2)
## 
## Call:
## lm(formula = y ~ x1)
## 
## 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 ***
## x1            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

Se observa que el proceso de selección de variables indica que el modelo que incluye solo diametro es suficiente para explicar el peso, ya que incorporar la altura no genera ningun beneficio adicional al modelo en terminos de ajuste.

Transformación de variables

Ya que el modelo estimado multiple indica que es suficiente solo usar el diametro, veamos que se puede hacer para mejorar el modelo solo con esta variable.

plot(x1,y)

plot(mod_multiple2)

Como se observa en las gráficas, logramos identificar que no existe una relación lineal entre el diametro con el peso del árbol, lo cual nos indica que podriamos mejorar el ajuste con una transformación. Se propone el logaritmo natural de y por su comportamiento creciente.

plot(x1, log(y))

mod_trans=lm(log(y)~x1)
summary(mod_trans)
## 
## Call:
## lm(formula = log(y) ~ x1)
## 
## 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 ***
## x1           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
plot(mod_trans)

Como se puede observar el modelo con la variable y transformada con logaritmo, mejora en ajuste a un 88.5% y adicionalmente los supuestos de linealidad y normalidad a nivel gráfico se ajustan mejor.

¿Que pasa al predecir con este modelo transformado, por ejemplo cual sería el peso estimado de un arbol con diametro 5.2m ?

y_log=predict(mod_trans, list(x1=5.2))
exp(y_log)
##        1 
## 16.03126

Un árbol con 5.2 m de diametro tiene un peso de 16 Ton.

Incluir variables predictoras de tipo categórico en el modelo

Vamos a incluir en el modelo la variable “finca”. Veamos como se interpretan estos resultados.

x3=datos_biomasa$finca
boxplot(y~x3,col="gray")

mod_trans_multi=lm(log(y)~x1+x3)
summary(mod_trans_multi)
## 
## Call:
## lm(formula = log(y) ~ x1 + x3)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.269702 -0.058652  0.000277  0.074438  0.233730 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.260860   0.045803  27.528  < 2e-16 ***
## x1          0.274018   0.009047  30.287  < 2e-16 ***
## x3FINCA_2   0.042002   0.031983   1.313    0.193    
## x3FINCA_3   0.227430   0.028554   7.965 6.23e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1069 on 86 degrees of freedom
## Multiple R-squared:   0.94,  Adjusted R-squared:  0.9379 
## F-statistic: 449.1 on 3 and 86 DF,  p-value: < 2.2e-16