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 acontinuación:

library(readxl)
datos_biomasa = read_excel("~/Desktop/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) y la altura como variable predictora (x). Se espere 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 por que el coeficiente de correlación de pearson es de 0.85, indicando que la altura del arbol 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 toneladas 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 variablidad del peso del arbol.

Con base en el valor p se observa que la variable altura es significativa en el modelo. Es decir la altura efectivamente es una 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 metros, si se estima que por cada tonelada de madera la compañia logra producir una cantidad de productos (papel, carton,…), cuyo valor estimado esta 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 una alto 1445 millones.

Modelo de Regresión Lineal Multiple

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

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_muliple1=lm(y~x1+x2)
summary(mod_muliple1)
## 
## 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 el diametro explica en parte lo de la altura ofrece al modelo (peso).

Selección de Variables

mod_multiple2=step(mod_muliple1)
## 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 graficas logramos identificar que no existe una relación lineal entre el diametro con el peso del arbol lo cual nos indica que podriamos mejorar el ajuste con una transformación. Se pronone 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 tranformada con logaritmo mejora en ajuste a un 88.5% y adicionalmente los supuestos de linealidad y normalidad a nivel grafico se ajustan mejor.

¿Que pasa al predecir con este modelo transformado, por ejemplo cual seria el peso estimado de un arbol con diametro de 5.2 metros?

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

Un arbol con 5.2 metros de diametro tiene un peso aproximado de 16 toneladas.

Incluir variables predictoras de tipo categorico en el modelo

Vamos 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

Observamos que la variable categorica finca es muy relevante en el modelo, ya que el ajuste del modelo se incrementa a 94% un valor muy alto. Adicionalmente se observan diferencias entre las fincas, en general la finca 1 presenta los arboles con el menor peso.

¿Cual seria el peso estimado de un arbol con diametro de 5.5 y que se encuentre en la finca 3?

y_trans2=predict(mod_trans_multi,list(x1=5.5,x3="FINCA_3"))
exp(y_trans2)
##        1 
## 19.99313
MAE2=mean(abs(exp(mod_trans_multi$residuals)))
MAE2
## [1] 1.00544

El modelo predice que el peso de ese arbol es de 19.99 aprox 20 toneladas. Adicionalmente podemos notar el MAE del modelo transformado es mucho menor que el modelo inicial solo con la altura.