Regresión simple y múltiple
Julieth Cerón
2023-05-10
Ejemplo modelo de árboles
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("D:/arboles.xlsx")
datos_biomasa## # A tibble: 90 × 5
## finca mg peso diametro altura
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 FINCA_1 GENOTIPO_1 13.7 4.7 5
## 2 FINCA_1 GENOTIPO_1 14.6 5.3 5.6
## 3 FINCA_1 GENOTIPO_1 15.9 4.8 5.8
## 4 FINCA_1 GENOTIPO_1 8.99 3.2 4.3
## 5 FINCA_1 GENOTIPO_1 6.99 2.2 3.3
## 6 FINCA_1 GENOTIPO_2 19.3 6.3 7.9
## 7 FINCA_1 GENOTIPO_2 21.4 6.6 8.3
## 8 FINCA_1 GENOTIPO_2 13.8 5.3 7.3
## 9 FINCA_1 GENOTIPO_2 11.9 4.9 6.7
## 10 FINCA_1 GENOTIPO_2 16.6 5.9 7.1
## # ℹ 80 more rowsPaso 1 - Plantear las variables del Modelo
Se desea evaluar la relación entre las variables peso 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
x=datos_biomasa$altura
plot(x,y)
cor(x,y)## [1] 0.8582009Se 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.891Como se observa el coeficiente B1 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-16Se 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 1444990De 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
x1=datos_biomasa$diametro
x2=datos_biomasa$altura
plot(datos_biomasa[,3:5])
cor(datos_biomasa[,3:5])## peso diametro altura
## peso 1.0000000 0.908123 0.8582009
## diametro 0.9081230 1.000000 0.9355360
## altura 0.8582009 0.935536 1.0000000mod_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-16Como 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 que 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.80summary(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-16Se 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)
Validación de supuestos
Podemos observar respecto a los supuestos sobre el error ei lo siguiente:
Media cero: Se cumple por defecto. Varianza Constante: Se observa en la grafica 1 de residuales vs ajustados que el comportamiento es aleatorio no con alguna tendencia en particular que indique problemas. Se valida grafiamente. Normalidad: Se observa en la grafica 2 que los datos se ajustan bien a la linea de normalidad en el qqplot. Es decir se valida graficamente. Independencia: Dado que estos registros no corresponden a datos en el tiempo no se tiene un orden temporal para realizar la validación de este supuesto. Se valida por definición del tipo de datos de corte transversal.
par(mfrow=c(2,2))
plot(mod_multiple2)
Como se observa en las graficas anteriores, 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 propone el logaritmo natural de y por su comportamiento creciente.
plot(x1,log(y))
Corremos el modelo con la nueva transformación:
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-16Validamos supuestos nuevamente:
par(mfrow=c(2,2))
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.03126Un 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")
Estimamos el modelo con la variable categórica:
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-16Observamos 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.
TAREA
1. ¿Cual seria el peso estimado de un arbol con diametro de 5.5 y que se encuentre en la finca 3?
2. ¿Cual seria el peso estimado de un arbol con diametro de 3.5 y que se encuentre en la finca 2?
3. Calcule el MAE del modelo transformado