Tarea regresion multivariada

en primer lugar se importan los datos del excel “arboles”

library(readxl)
datos_biomasa = read_excel("C:/Users/migue/Desktop/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 rows

se definen las variables independientes: diametro, altura y finca
se define como variable dependiente el peso

se grafican las relaciones entre las variables cuantitativas, haciendo la diferenciacion teniendo en cuenta la variable finca

y = datos_biomasa$peso
x1 = datos_biomasa$diametro
x2 = datos_biomasa$altura
x3 = datos_biomasa$finca

#finca 1
plot(datos_biomasa[1:30,3:5], main = "FINCA 1")

#finca 2
plot(datos_biomasa[31:60,3:5], main = "FINCA 2")

#finca 3
plot(datos_biomasa[61:90,3:5], main = "FINCA 3")

#diferencia entre fincas
boxplot(y~x3, col = "gray")

el boxplot muestra las notables diferencias que hay entre fincas, por lo que esta variable categorica sera importante en la prediccion mas adelante

se muestran las correlaciones entre las variables, nuevamente teniendo en cuenta la separacion segun la finca

#finca 1
cor(datos_biomasa[1:30,3:5])
##               peso  diametro    altura
## peso     1.0000000 0.9313303 0.8198852
## diametro 0.9313303 1.0000000 0.9441192
## altura   0.8198852 0.9441192 1.0000000
#finca 2
cor(datos_biomasa[31:60,3:5])
##               peso  diametro    altura
## peso     1.0000000 0.9337049 0.8754702
## diametro 0.9337049 1.0000000 0.9154492
## altura   0.8754702 0.9154492 1.0000000
#finca 3
cor(datos_biomasa[61:90,3:5])
##               peso  diametro    altura
## peso     1.0000000 0.9082605 0.9082983
## diametro 0.9082605 1.0000000 0.9480352
## altura   0.9082983 0.9480352 1.0000000

segun el ejercicio hecho en clase se encontro que la variable diametro (x1) es la mas adecuada para explicar el peso de los arboles.

Se muestra a continuacion, la grafica entre diametro y peso para cada finca

plot(x1,y, main = "FINCAS 1, 2 y 3")

plot(x1[1:30],y[1:30], main = "FINCA 1")

plot(x1[31:60],y[31:60], main = "FINCA 2")

plot(x1[61:90],y[61:90], main = "FINCA 3")

se crea el modelo de regresion simple con la variable diametro para tener una base de comparacion mas adelante con los distintos modelos

mod_simple=lm(y~x1)
summary(mod_simple)
## 
## 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

la variable diametro explica es un 82.3% la variable peso

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

se grafica para validar los supuestos y se observa que la varianza no es constante, ademas viendo las graficas de diametro vs peso, se observa una tendencia un poco curva, por lo tanto se plantea una trasnformacion de logaritmo natural del peso

plot(x1,log(y), main = "FINCAS 1, 2 y 3")

plot(x1[1:30],log(y[1:30]), main = "FINCA 1")

plot(x1[31:60],log(y[31:60]), main = "FINCA 2")

plot(x1[61:90],log(y[61:90]), main = "FINCA 3")

se grafica con la transformada y ahora si se aprecia mejor la linealidad tanto en el grupo completo de datos, como en la segmentacion por fincas

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

se crea el modelo con la transformada y se observa que explica un 88.5% el peso, mejor que en el caso anterior

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

tambien se verifican los supuestos y se nota que la varianza se comporta mas estable a lo largo de todos los datos

finalmente se crea el modelo con las variables diametro y finca

mod_multiple = lm(log(y)~x1+x3)
summary(mod_multiple)
## 
## 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

se confirma una mejora en el valor R squared, este modelo multivariado explica un 93.8% la variable peso

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

se verifican los supuestos y se observa una varianza constante y normalidad con el qqplot

finalmente se muestran las predicciones solicitadas en la tarea

y_pred1 = predict(mod_multiple,list(x1 = 5.5, x3 = "FINCA_3"))
exp(y_pred1)
##        1 
## 19.99313
y_pred2 = predict(mod_multiple,list(x1 = 3.5, x3 = "FINCA_2"))
exp(y_pred2)
##        1 
## 9.601551

el peso estimado de un arbol con diametro 5.5 en la finca 3 es de: 19.99 tn

el peso estimado de un arbol con diametro 3.5 en la finca 2 es de: 9.60 tn

finalmente se muestra el MAE tanto del modelo con transformada como el modelo multivariado, siendo este ultimo el que tiene un MAE menor

#MAE modelo transformado
MAEt = mean(abs(mod_trans$residuals))
MAEt
## [1] 0.1179231
#MAE modelo multivariado
MAEm = mean(abs(mod_multiple$residuals))
MAEm
## [1] 0.08240871