Informe 1 – Modelo de Regresión Lineal Simple y Múltiple
Carlos Luis Mora Cañas
2022-08-21
Taller Ciencia de Datos
Analizar el comportamiento de los precios de las Acciones de Ecopetrol según la variación del precio del barril de petróleo WTI producido en Colombia. Se tienen los siguientes precios.
acciones <- read_csv("C:/Users/icm2363a/Documents/R/Informe-1/datasets/acciones.csv",show_col_types = FALSE)
head(acciones)## # A tibble: 6 × 3
## Fecha Accion_Ecopetrol_COP Petroleo_WTI_Dolares
## <chr> <dbl> <dbl>
## 1 14/12/2015 1090 35.6
## 2 15/12/2015 1170 36.3
## 3 16/12/2015 1160 37.4
## 4 17/12/2015 1230 35.0
## 5 21/12/2015 1155 34.5
## 6 22/12/2015 1165 35.8A. Proponga un modelo de regresión lineal simple que permita predecir el valor de las Acciones de Ecopetrol con base en el Precio del barril de petróleo en Colombia. Indique la ecuación de regresión y el valor del R2
plot(acciones$Accion_Ecopetrol_COP,acciones$Petroleo_WTI_Dolares)
#### Estimacion Modelo MCO
modelo_acciones <-(lm(formula = (acciones$Accion_Ecopetrol_COP ~ acciones$Petroleo_WTI_Dolares), data = acciones))
summary(modelo_acciones)##
## Call:
## lm(formula = (acciones$Accion_Ecopetrol_COP ~ acciones$Petroleo_WTI_Dolares),
## data = acciones)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.90 -40.74 -15.94 33.40 136.82
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 177.768 232.828 0.764 0.45627
## acciones$Petroleo_WTI_Dolares 26.192 6.542 4.004 0.00102 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.13 on 16 degrees of freedom
## Multiple R-squared: 0.5005, Adjusted R-squared: 0.4692
## F-statistic: 16.03 on 1 and 16 DF, p-value: 0.001024La ecuación de la regresión lineal es y^ = 177.768 + 26.192x^
El valor de R2 es 0.50
================================================================================
B. Pruebe la significancia del modelo propuesto en “a)” plantee las hipótesis respectivas y use el concepto de Valor _p para tomar la decisión sobre las hipótesis. Use α = 0.05
La hipotesis nula donde b1 es igual 0 y que b1 es diferente de 0. El p-Value es equivalente a 0.001024 que es menor a α = 0.05 y podemos rechazar la hipotesis nula o donde b1 es diferente de 0.
================================================================================
C. Haga un análisis de los residuos. ¿Qué supuesto no se cumple?
Para todos los casos se toma nivel de significancia 0.05.
Distribución normal
# Residuales
residuales_acciones <- modelo_acciones$residuals
shapiro.test(residuales_acciones)
plot2 <- ggplot(data = acciones, aes(Accion_Ecopetrol_COP, modelo_acciones$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = acciones, aes(Petroleo_WTI_Dolares, modelo_acciones$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot2, plot3)
No se cumple el supuesto de normalidad, dado que el p-Value es 0.04276, dando a decir que hay significancia para rechaza dicha hipotesis nula equivalente a normal.
Varianza constante - homocedasticidad
lmtest::bptest(modelo_acciones)##
## studentized Breusch-Pagan test
##
## data: modelo_acciones
## BP = 0.029563, df = 1, p-value = 0.8635bptest(modelo_acciones) ##
## studentized Breusch-Pagan test
##
## data: modelo_acciones
## BP = 0.029563, df = 1, p-value = 0.8635El supuesto de homocedasticidad da un p-value es 0.86, lo que nos indica que no hay significancia para rechaza la hipotesis nula equivalente a varianza constante.
Errores independientes
lmtest::dwtest(modelo_acciones) ##
## Durbin-Watson test
##
## data: modelo_acciones
## DW = 0.74504, p-value = 0.0004666
## alternative hypothesis: true autocorrelation is greater than 0El supuesto de los errores independientes da un p-value de 0.0004, lo que nos indica que hay evidencia para rechaza la hipotesis nula equivalente a no autocorrelación.
================================================================================
E. Concluya sobre la validez del modelo propuesto en a)
Dicho modelo no cumple 2 supuestos de la regresión lineal, por este caso que los errores sean normales y que no exista autocorrelación, ademas el R2 no se ajusta a los datos, todo lo anterior mencionado hace que en el modelo no se pueda establecer conclusiones certeras; Se recomienda hacer transformaciones a los datos para que se parezcan más a un modelo lineal y tratar de obtener más muestras
================================================================================
Valores SSMLV - Inflación
Los siguientes datos corresponden a la INFLACION y al SALARIO MINIMO LEGAL MENSUAL (SMLM) desde el año 1999 para Colombia.
salario <- read_csv("C:/Users/icm2363a/Documents/R/Informe-1/datasets/smmlv.csv",show_col_types = FALSE)
head(salario,3)## # A tibble: 3 × 3
## AÑO INFLACION SMLM
## <chr> <dbl> <dbl>
## 1 1999 9.23 236460
## 2 2000 8.75 260100
## 3 2001 7.65 286000A. Escriba la ecuación del modelo de regresión lineal simple
# estimacion del modelo por MCO
modelo_salario <- lm(salario$SMLM ~ salario$INFLACION)
summary(modelo_salario)##
## Call:
## lm(formula = salario$SMLM ~ salario$INFLACION)
##
## Residuals:
## Min 1Q Median 3Q Max
## -758977 -509891 -363360 -93378 5585918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 478201 943505 0.507 0.620
## salario$INFLACION 56039 162484 0.345 0.735
##
## Residual standard error: 1507000 on 15 degrees of freedom
## Multiple R-squared: 0.007867, Adjusted R-squared: -0.05827
## F-statistic: 0.1189 on 1 and 15 DF, p-value: 0.735La ecuación de la regresión lineal es y^ = 478201 + 56039x^
================================================================================
B. Plantee y valide las hipótesis correspondientes a la linealidad general del modelo propuesto en a)
plot(salario$SMLM, salario$INFLACION)
Nivel de significancia 0.05 que es el usual.
Distribución normal
residuales_salario <- modelo_salario$residuals
shapiro.test(residuales_salario)##
## Shapiro-Wilk normality test
##
## data: residuales_salario
## W = 0.41241, p-value = 2.546e-07No se cumple con el supuesto de normalidad dado al p-value es 0.000007, lo cual indica que hay significancia para rechaza la hipotesis nula equivalente a normal.
Varianza constante
lmtest::bptest(modelo_salario) ##
## studentized Breusch-Pagan test
##
## data: modelo_salario
## BP = 0.54706, df = 1, p-value = 0.4595bptest(modelo_salario)##
## studentized Breusch-Pagan test
##
## data: modelo_salario
## BP = 0.54706, df = 1, p-value = 0.4595Cumplimos el supuesto de homocedasticidad porque nuestro valor p es 0.45, indica que no hay significancia para rechaza la hipotesis nula equivalente a varianza constante.
Errores independientes
lmtest::dwtest(modelo_salario) ##
## Durbin-Watson test
##
## data: modelo_salario
## DW = 0.94387, p-value = 0.003993
## alternative hypothesis: true autocorrelation is greater than 0dwtest(modelo_salario)##
## Durbin-Watson test
##
## data: modelo_salario
## DW = 0.94387, p-value = 0.003993
## alternative hypothesis: true autocorrelation is greater than 0No se cumple con el supuesto de errores independientes, porque el p-value es equivalente a 0.003, dando a decir que hay evidencia para rechaza la hipotesis nula equivalente a no autocorrelación.
================================================================================
C. Indique e interprete el coeficiente de correlación del modelo propuesto en a)
EL coeficiente de correlación es equivalente a la raiz cuadrada del coeficiente de determinación R2
sqrt(0.007867)## [1] 0.08869611Se Observa poco correlacionados (8.8%)
================================================================================
D. Interprete cada uno de los coeficientes del modelo propuesto en a)
El coeficiente de determinación indica que la pendiente del modelo se acomoda poco a los datos debido a que el r2 es equvialente a 0.007867.
El p-value equivalente a 0.735 no evidencia que el modelo sea significativo.
El intercepto estimado evidencia que los valores comienzan en 478.201 COP, pero que su p-value no rechaza la hipotesis nula.
El coeficiente de inflación estimado evidencia que su coeficiente es 56.039, pero que su p-value no rechaza la hipotesis nula, el modelo no puede explicarse por la variable inflación.
================================================================================
E. Construya una gráfica de residuales y haga un análisis cualitativo de los supuestos del modelo propuesto en a)
par(mfrow=c(2,2))
plot(modelo_salario)
qqnorm(modelo_salario$residuals)
qqline(modelo_salario$residuals)
Los graficos residuals vs fitted se observa que no existe una aleatoriedad de los errores, es decir, no son indepedinetes. Las escalas que se estan comparando no permiten ver con claridad si existe normalidad porque su recta no esta tan marcada, en cuanto a la homocedasticidad evidenciamos que la recta no es horizontal sino tiene una relación directa entre la escala y la ubicación. Se recomienda hacer primero una transformación para realizar una interpretación más adecuada de los gráficos.
================================================================================
F. Comente sobre la conveniencia de usar el modelo propuesto en a) para predecir el SMLM para Colombia
No es conveniente realizar predicciones sobre el modelo, se indica que el modelo con la variable independiente Inflación no puede explicar los cambios en el salario minimo.
================================================================================
Precios vivienda
Con base en los datos de precios de vivienda de la actividad en clase realizar un informe que contenga los siguientes puntos utilizando R y RMarkdown.
A. Realice un filtro a la base de datos e incluya solo las ofertas de apartamentos, de la zona norte de la ciudad con precios inferiores a los 500 millones de pesos y áreas menores a 300 mt2. Presente los primeros 3 registros de la base y algunas tablas que comprueben la consulta. (Adicional un mapa con los puntos de la base, discutir si todos los puntos se ubican en la zona norte o se presentan valores en otras zonas, por que?)
viviendas <- read_excel("C:/Users/icm2363a/Documents/R/Informe-1/datasets/datos_vivienda.xlsx")
filtro1 <- subset(viviendas, viviendas$Tipo == "Apartamento" & viviendas$precio_millon < 500 &
viviendas$Zona == "Zona Norte" & viviendas$Area_contruida < 300)
head(filtro1, 3)## # A tibble: 3 × 12
## Zona piso Estrato preci…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio corde…⁵
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <chr> <chr> <dbl>
## 1 Zona… 2 3 135 56 1 1 3 Apar… torre… -76.5
## 2 Zona… NA 3 78 54 2 1 3 Apar… chimi… -76.5
## 3 Zona… NA 5 340 106 2 2 3 Apar… la fl… -76.5
## # … with 1 more variable: Cordenada_latitud <dbl>, and abbreviated variable
## # names ¹precio_millon, ²Area_contruida, ³parqueaderos, ⁴Habitaciones,
## # ⁵cordenada_longitud#Verificación
unique(filtro1$Zona)## [1] "Zona Norte"unique(filtro1$Tipo)## [1] "Apartamento"leaflet(filtro1) %>% addTiles() %>% addProviderTiles(providers$CartoDB.Positron)%>%
addCircleMarkers(
lng = filtro1$cordenada_longitud,
lat = filtro1$Cordenada_latitud,
label=paste(filtro1$Barrio,filtro1$precio_millon,sep=":"),
radius= filtro1$precio_millon/100)La mayoría de puntos se encuentran en zona norte, sin embargo, si se evidencia puntos distribuidos en la zona centro y sur, se resalta que algunos marcan que el barrio es Calí y Acopi, en cuyo caso son barrios donde la zona y el barrio n ser ajustados. Se resalta también un cluster muy importante en el Barrio Lili
================================================================================
B. Realice un análisis exploratorio de datos enfocado en la correlación entre la variable respuesta (precio del apartamento) en función del área construida, estrato y si tiene parqueadero. Use gráficos interactivos con plotly e interprete los resultados.
g1 <- ggplot(data =filtro1, aes(y=Area_contruida, x=precio_millon, colour=parqueaderos)) + geom_point() + labs(title ="Precio por área construida", x="Precio Millón", y="Mts2")
ggplotly(g1)Se evidencia una correlación directa entre el precio de los apartamentos con el área construida, aunque en un ánalisis exploratorio inicial podemos ver que su varianza no es constante. Otro punto importante sobre la presente gráfica es que también se evidencia relación directa del precio con el número de apartamento, el precio para un apartamento sin parqueaderos se encuentra en su mayoría en menos de 100 millones, mientras que cuando tienen 2 parqueaderos en su mayoría no bajan de 250 millones
g2 <- ggplot(data =filtro1, aes(y=Estrato, x=precio_millon))+ geom_smooth(alpha=0.3) + geom_point(alpha=0.6) + labs(title ="Precio por estrato", x="Precio Millón", y="Estrato")
ggplotly(g2)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'En menor medida se nota la relación directa que existe entre el estrato y el precio, sin embargo, logramos ver que a medida que aumenta el estrato si hay rangos superiores en los datos, mientras la mayoría de datos para el estrato 3 ronda los precios de 65 a 200 millones, el estrato 5 tiene precios que ronda 130 a 500 millones.
================================================================================
C. Estime un modelo de regresión lineal múltiple con las variables del punto anterior e interprete los coeficientes si son estadísticamente significativos. Las interpretaciones deber están contextualizadas y discutir si los resultados son lógicos. Adicionalmente interprete el coeficiente R2 y discuta el ajuste del modelo e implicaciones (que podrían hacer para mejorarlo).
modelo <- lm(precio_millon ~ Area_contruida + parqueaderos + Estrato, data= filtro1)
summary(modelo)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + parqueaderos +
## Estrato, data = filtro1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -199.038 -31.160 -0.589 28.955 231.217
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -137.72572 8.68837 -15.852 < 2e-16 ***
## Area_contruida 0.89025 0.06182 14.402 < 2e-16 ***
## parqueaderos2 41.29670 4.96122 8.324 2.57e-16 ***
## parqueaderos3 42.45052 27.06700 1.568 0.11710
## parqueaderos4 79.15804 53.42927 1.482 0.13875
## parqueaderosNA -10.70156 3.88228 -2.757 0.00594 **
## Estrato 69.27921 2.25772 30.685 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.32 on 1070 degrees of freedom
## Multiple R-squared: 0.7673, Adjusted R-squared: 0.766
## F-statistic: 588.2 on 6 and 1070 DF, p-value: < 2.2e-16Para comenzar notamos que la variable área construida tiene significancia representativa para explicar el precio, así mismo el estrato tiene un p valor significativo para explicar el modelo, sin embargo, de los parqueaderos podemos ver que solo sin parqueadero y con 2 parqueadero logran explicar la variable dependiente, esto es debido a que para el caso de 1 parqueadero existen puntos a lo largo de los diferentes valores y no se conforman clusters, para el caso de 3 y 4 parqueaderos se obtienen muy pocos valores. Se recomienda que al igual que los parqueaderos miremos los efectos de dividir en dummis los estratos porque también en el análisis exploratorio se evidencio un amplio rango para el caso del estrato 5.
Respecto al coeficiente de determinación notamos un ajuste de los datos aceptable, sin embargo, se recomienda ajustar aún más para ello se pueden proponer transformación de los datos que permitan parecer más los datos a una recta, acercarse a una varianza constante y revisar si existe una mejora en el R2, por otro lado como lo vimos en el mapa, se evidencian zonas diferentes a la del norte, que incluye puntos de zona sur y zona centro, por tal razón se recomienda también hacer un filtro de estos puntos para que los valores restantes sean geograficamente homogeneos
================================================================================
D. Realice la validación de supuestos del modelo e interprete los resultados (no es necesario corregir en caso de presentar problemas solo realizar sugerencias de que se podría hacer).
#Distribución normal
residuales_viviendas <- modelo$residuals
shapiro.test(residuales_viviendas)##
## Shapiro-Wilk normality test
##
## data: residuales_viviendas
## W = 0.98724, p-value = 4.471e-08No cumplimos con el supuesto de normalidad porque nuestro valor p es 4.471e-08, indica que hay significancia para rechaza la hipotesis nula equivalente a normal.
#Varianza constante
lmtest::bptest(modelo) ##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 154.83, df = 6, p-value < 2.2e-16No cumplimos el supuesto de homocedasticidad porque nuestro valor p es 2.2e-16, indica que hay significancia para rechaza la hipotesis nula equivalente a varianza constante, esto es validado por la gráfica inicial
#Errores independientes
lmtest::dwtest(modelo) ##
## Durbin-Watson test
##
## data: modelo
## DW = 1.7497, p-value = 1.766e-05
## alternative hypothesis: true autocorrelation is greater than 0No cumplimos con el supuesto de errores independientes porque nuestro valor p es equivalente a 0.003, indica que hay evidencia para rechaza la hipotesis nula equivalente a no autocorrelación.
Sin duda, no cumplimos ninguno de los supuestos, por lo tanto es limitada la inferencia que se puede hacer de los coeficientes del modelo, sin embargo, luego del análisis exploratorio y los resultados de la regresión evidenciamos que se podría encontrar la forma de la pendiente.
================================================================================ #### E. Con el modelo identificado predecir el precio de un apartamento con 100 mt2, de estrato 4 y con parqueadero. ¿Si este apartamento lo están ofreciendo en 450 millones cual seria su opinión con base en el resultado del modelo considera que es una buena oferta?
modelo##
## Call:
## lm(formula = precio_millon ~ Area_contruida + parqueaderos +
## Estrato, data = filtro1)
##
## Coefficients:
## (Intercept) Area_contruida parqueaderos2 parqueaderos3 parqueaderos4
## -137.7257 0.8903 41.2967 42.4505 79.1580
## parqueaderosNA Estrato
## -10.7016 69.2792input <- data.frame(Area_contruida=c(100), parqueaderos=c("1"), Estrato=c(4))
predict(modelo,
newdata=input,
interval="confidence",
level=0.95) ## fit lwr upr
## 1 228.4162 223.1704 233.662Si el precio de ese apartamento es de 450 millones, estan cobrando el doble de lo que vale en el mercado, no tomaría esa opción de compra y me acogería a alternativas, mala oferta.
================================================================================ #### F. Con las predicciones del modelo sugiera potenciales ofertas para una persona interesada en un apartamento en la zona norte con mas de 100 mt2 de área, de estrato 4, que tenga parqueadero y tenga encuentra que la persona tiene un crédito preaprobado de máximo 400 millones de pesos. Realice un análisis y presente en un mapa al menos 5 ofertas potenciales que debe discutir.
input <- data.frame(Area_contruida=c(265), parqueaderos=c("1"), Estrato=c(4))
predict(modelo,
newdata=input,
interval="confidence",
level=0.95) ## fit lwr upr
## 1 375.3076 352.1567 398.4585Dada la situación del mercado, el cliente podría obtener un apartamento con esas características que ronde entre los 100 y 265 metros cuadrados. A continuación presentamos las ofertas más relevantes de acuerdo a los parametros solicitados, se recomienda ver los circulos más grandes quienes contienen las áreas más significativas.
filtro2 <- subset(viviendas, viviendas$Tipo == "Apartamento" & viviendas$precio_millon < 400 &
viviendas$Zona == "Zona Norte" & viviendas$Area_contruida < 265 &
viviendas$Area_contruida > 100 & viviendas$Estrato == 4)
leaflet(filtro2) %>% addTiles() %>% addProviderTiles(providers$CartoDB.Positron)%>%
addCircleMarkers(
lng = filtro2$cordenada_longitud,
lat = filtro2$Cordenada_latitud,
label=paste('Barrio: ',filtro2$Barrio,'Precio: ',filtro2$precio_millon, 'Estrato: ',
filtro2$Estrato,'Metros: ',filtro2$Area_contruida, sep=' '),
radius= filtro2$Area_contruida/50)================================================================================
Datos arboles
Con base en los datos de arboles proponga un modelo de regresión lineal múltiple que permita predecir el peso del árbol en función de las covariables que considere importantes y seleccionándolas de acuerdo con un proceso adecuado. Tenga en cuenta realizar una evaluación de la significancia de los parámetros, interpretación y proponga un método de evaluación por medio de validación cruzada. Presente métricas apropiadas como el RMSE y MAE.
arboles <- read_excel("C:/Users/icm2363a/Documents/R/Informe-1/datasets/data arboles.xlsx")
head(arboles,3)## # A tibble: 3 × 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.8Descriptiva de los datos
dataf.2= arboles
str(dataf.2)## tibble [90 × 5] (S3: tbl_df/tbl/data.frame)
## $ finca : chr [1:90] "FINCA_1" "FINCA_1" "FINCA_1" "FINCA_1" ...
## $ mg : chr [1:90] "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" ...
## $ peso : num [1:90] 13.73 14.58 15.88 8.99 6.99 ...
## $ diametro: num [1:90] 4.7 5.3 4.8 3.2 2.2 6.3 6.6 5.3 4.9 5.9 ...
## $ altura : num [1:90] 5 5.6 5.8 4.3 3.3 7.9 8.3 7.3 6.7 7.1 ...summary(dataf.2)## finca mg peso diametro
## Length:90 Length:90 Min. : 5.98 Min. :2.200
## Class :character Class :character 1st Qu.:13.64 1st Qu.:4.525
## Mode :character Mode :character Median :17.48 Median :5.400
## Mean :18.77 Mean :5.446
## 3rd Qu.:22.80 3rd Qu.:6.500
## Max. :47.87 Max. :8.800
## altura
## Min. : 3.300
## 1st Qu.: 5.225
## Median : 6.450
## Mean : 6.634
## 3rd Qu.: 7.875
## Max. :11.300describe(dataf.2)## vars n mean sd median trimmed mad min max range skew kurtosis
## finca* 1 90 2.00 0.82 2.00 2.00 1.48 1.00 3.00 2.00 0.00 -1.53
## mg* 2 90 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.00 0.00 -2.02
## peso 3 90 18.77 8.16 17.48 17.96 6.35 5.98 47.87 41.89 1.11 1.55
## diametro 4 90 5.45 1.45 5.40 5.46 1.48 2.20 8.80 6.60 -0.07 -0.36
## altura 5 90 6.63 1.80 6.45 6.55 1.93 3.30 11.30 8.00 0.38 -0.44
## se
## finca* 0.09
## mg* 0.05
## peso 0.86
## diametro 0.15
## altura 0.19# Diagrama de Dispersion
g1 <- ggplot(data =arboles, aes(y=diametro, x=peso, colour=finca)) + geom_point() + labs(title ="Peso por diametro", x="Peso", y="Diametro")
ggplotly(g1)# Histogramas
histograma <- ggplot(dataf.2, aes(x=peso)) +
ggtitle("peso") +
geom_histogram(color="#28324a", fill="#3c78d8")
histograma## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
histograma1 <- ggplot(dataf.2, aes(x=diametro)) +
ggtitle("diametro") +
geom_histogram(color="#28324a", fill="#3c78d8")
histograma1## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
histograma2 <- ggplot(dataf.2, aes(x=altura)) +
ggtitle("altura") +
geom_histogram(color="#28324a", fill="#3c78d8")
histograma2## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#### graficas descriptivas
ggpairs(dataf.2, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Observamos una relación lineal directa muy marcada, incluso para los 3
tipos de fincas e inicialmente con una varianza constante, entre más
aumenta el peso, más alto es el diametro.
g2 <- ggplot(data =arboles, aes(y=altura, x=peso, colour=mg)) + geom_point() + labs(title ="Peso por altura", x="Peso", y="Altura")
ggplotly(g2)#Diagrama de dispersion
plot(arboles$peso, arboles$diametro,main = "Diagrama de disperión",xlab = "Peso", ylab = "Diametro")
plot(arboles$peso, arboles$altura,main = "Diagrama de disperión",xlab = "Peso", ylab = "Altura")
pairs(arboles$peso ~ arboles$diametro)
pairs(arboles$peso ~ arboles$altura)
De igual manera, la altura presenta una relación directa marcada, los puntos también parecen estar un canal con varianza constante tanto para genotipo 1 y 2.
No parece tener multicolinealidad las 2 variables, dado que sus niveles de correlación varian respecto al peso
Regresión inicial -> Análisis de significancia de las variables, supuestos y redimiento inicial del modelo
modelo_arboles1 <- lm(peso ~ finca + mg + diametro + altura, data=arboles)
summary(modelo_arboles1)##
## Call:
## lm(formula = peso ~ finca + mg + diametro + altura, data = arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1009 -1.8569 -0.5094 1.5578 12.8691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -13.95177 1.68295 -8.290 1.59e-12 ***
## fincaFINCA_2 -0.03095 0.99140 -0.031 0.975166
## fincaFINCA_3 3.51938 0.83466 4.217 6.23e-05 ***
## mgGENOTIPO_2 -4.50270 1.23667 -3.641 0.000468 ***
## diametro 2.57058 0.76282 3.370 0.001138 **
## altura 2.98566 0.76616 3.897 0.000195 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.983 on 84 degrees of freedom
## Multiple R-squared: 0.8738, Adjusted R-squared: 0.8662
## F-statistic: 116.3 on 5 and 84 DF, p-value: < 2.2e-16modelo <- (lm(formula = (peso ~ altura + diametro), data = arboles))
summary(modelo)##
## Call:
## lm(formula = (peso ~ altura + diametro), data = arboles)
##
## 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 ***
## altura 0.3132 0.5751 0.544 0.587
## diametro 4.7395 0.7128 6.649 2.49e-09 ***
## ---
## 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#Distribución normal
residuales_arboles <- modelo_arboles1$residuals
shapiro.test(residuales_arboles)##
## Shapiro-Wilk normality test
##
## data: residuales_arboles
## W = 0.90738, p-value = 8.795e-06plot2 <- ggplot(data = arboles, aes(altura, modelo_arboles1$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = arboles, aes(diametro, modelo_arboles1$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot2, plot3)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Distribución normal de los residuos:
qqnorm(modelo_arboles1$residuals)
qqline(modelo_arboles1$residuals)
#No cumplimos con el supuesto de normalidad porque nuestro valor p es 8.795e-06, indica que hay significancia para rechaza la hipotesis nula equivalente a normal.
#Varianza constante
lmtest::bptest(modelo_arboles1) ##
## studentized Breusch-Pagan test
##
## data: modelo_arboles1
## BP = 10.437, df = 5, p-value = 0.06375#Cumplimos el supuesto de homocedasticidad porque nuestro valor p es 0.06375, indica que hay noy significancia para rechaza la hipotesis nula equivalente a varianza constante, esto es validado por la gráfica inicial donde no se nota una apertura del canal de los puntos.
# Errores independientes
lmtest::dwtest(modelo_arboles1) ##
## Durbin-Watson test
##
## data: modelo_arboles1
## DW = 1.4259, p-value = 0.0008882
## alternative hypothesis: true autocorrelation is greater than 0# Grafico
ggplot(data = arboles, aes(modelo_arboles1$fitted.values, modelo_arboles1$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#varianza constante homocedasticidad
bptest(modelo_arboles1) ##
## studentized Breusch-Pagan test
##
## data: modelo_arboles1
## BP = 10.437, df = 5, p-value = 0.06375# Matriz de correlación entre predictores.corrplot(cor(dplyr::select(arboles, altura, diametro)),
method = "number", tl.col = "black")
# Análisis de Inflación de Varianza (VIF):
vif(modelo_arboles1)## GVIF Df GVIF^(1/(2*Df))
## finca 1.874507 2 1.170097
## mg 3.866072 1 1.966233
## diametro 12.263683 1 3.501954
## altura 19.004481 1 4.359413# Autocorrelación
dwt(modelo_arboles1, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.2836953 1.425949 0.004
## Alternative hypothesis: rho != 0# Identificación de posibles valores atípicos o influyentes
outlierTest(modelo_arboles1)## rstudent unadjusted p-value Bonferroni p
## 53 5.243166 1.1819e-06 0.00010637
## 55 3.874839 2.1230e-04 0.01910700# Validar si es influyente
summary(influence.measures(modelo_arboles1))## Potentially influential observations of
## lm(formula = peso ~ finca + mg + diametro + altura, data = arboles) :
##
## dfb.1_ dfb.fFINCA_2 dfb.fFINCA_3 dfb.mGEN dfb.dmtr dfb.altr dffit cov.r
## 53 -0.25 0.15 -0.52 0.64 1.37_* -1.06_* 1.84_* 0.22_*
## 55 -1.15_* -0.19 -0.06 -0.80 -0.50 0.93 1.48_* 0.45_*
## 58 0.06 0.05 -0.06 0.07 0.20 -0.19 0.24 1.24_*
## cook.d hat
## 53 0.43 0.11
## 55 0.31 0.13
## 58 0.01 0.15influencePlot(modelo_arboles1)
## StudRes Hat CookD
## 34 -1.296678 0.1578091 0.05208659
## 40 1.133418 0.1640228 0.04186676
## 53 5.243166 0.1097032 0.42921387
## 55 3.874839 0.1274908 0.31336782No cumplimos con el supuesto de errores independientes porque nuestro valor p es equivalente a 0.003, indica que hay evidencia para rechaza la hipotesis nula equivalente a no autocorrelación.
# Metricas iniciales de rendimiento
# division de dataset
split = sample.split(arboles$peso, SplitRatio = 0.8)
nltrain = subset(arboles, split == TRUE)
nltest = subset(arboles, split == FALSE)
division <- arboles$peso %>% createDataPartition (p=0.8, list=FALSE)
entrena <- arboles [division,]
prueba <- arboles [-division,]
# Validaciones dataset
dim(entrena)## [1] 74 5dim(prueba)## [1] 16 5str(entrena)## tibble [74 × 5] (S3: tbl_df/tbl/data.frame)
## $ finca : chr [1:74] "FINCA_1" "FINCA_1" "FINCA_1" "FINCA_1" ...
## $ mg : chr [1:74] "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" ...
## $ peso : num [1:74] 13.73 14.58 15.88 8.99 19.34 ...
## $ diametro: num [1:74] 4.7 5.3 4.8 3.2 6.3 6.6 4.9 5.9 2.5 4.7 ...
## $ altura : num [1:74] 5 5.6 5.8 4.3 7.9 8.3 6.7 7.1 3.7 5.7 ...str(prueba)## tibble [16 × 5] (S3: tbl_df/tbl/data.frame)
## $ finca : chr [1:16] "FINCA_1" "FINCA_1" "FINCA_1" "FINCA_1" ...
## $ mg : chr [1:16] "GENOTIPO_1" "GENOTIPO_2" "GENOTIPO_1" "GENOTIPO_2" ...
## $ peso : num [1:16] 6.99 13.81 7.47 18.69 10.33 ...
## $ diametro: num [1:16] 2.2 5.3 2.2 6.3 3.6 5.6 4.3 6.6 6.4 7.7 ...
## $ altura : num [1:16] 3.3 7.3 3.5 8.1 4.6 6.2 6.3 8.6 7.4 10.1 ...##
dim(nltrain)## [1] 72 5dim(nltest)## [1] 18 5str(nltrain)## tibble [72 × 5] (S3: tbl_df/tbl/data.frame)
## $ finca : chr [1:72] "FINCA_1" "FINCA_1" "FINCA_1" "FINCA_1" ...
## $ mg : chr [1:72] "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" ...
## $ peso : num [1:72] 13.73 14.58 15.88 8.99 6.99 ...
## $ diametro: num [1:72] 4.7 5.3 4.8 3.2 2.2 6.3 6.6 4.9 5.9 3.1 ...
## $ altura : num [1:72] 5 5.6 5.8 4.3 3.3 7.9 8.3 6.7 7.1 4 ...str(nltest)## tibble [18 × 5] (S3: tbl_df/tbl/data.frame)
## $ finca : chr [1:18] "FINCA_1" "FINCA_1" "FINCA_1" "FINCA_1" ...
## $ mg : chr [1:18] "GENOTIPO_2" "GENOTIPO_1" "GENOTIPO_1" "GENOTIPO_1" ...
## $ peso : num [1:18] 13.81 7.03 15.83 7.47 12 ...
## $ diametro: num [1:18] 5.3 2.5 4.7 2.2 4.9 4.2 3 3.7 4.8 4.7 ...
## $ altura : num [1:18] 7.3 3.7 5.7 3.5 7 5.2 4.8 4.4 5.5 5.8 ...Creación del modelo
modelo_arboles <- lm(peso ~ finca + mg + diametro + altura, data=arboles)
#gráfica de regresiones en parejas, con línea de regresión variables numericas
datos1 <- arboles[,-1] # Elimina la primera columna
datos2 <- datos1[,-1] # Elimina la primera columna
scatterplotMatrix(datos2, peso ~ altura + diametro,
smooth = list(lty = 2), id = TRUE,
regLine = list(lty = 1, col = "red"),
col = "blue")## Warning in applyDefaults(legend, defaults = list(coords = NULL, pt.cex = cex, :
## unnamed legend arguments, will be ignored
aic <- AIC(modelo)
sprintf("AIC = %.2f", aic)## [1] "AIC = 483.20"#Estimación
estimaciones <- modelo_arboles %>% predict(prueba)
# prueba de homocedasticidad
ncvTest(modelo_arboles)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 23.29964, Df = 1, p = 1.3863e-06spreadLevelPlot(modelo_arboles)
##
## Suggested power transformation: 1.09914#guardar resultados
resultados <- data.frame (Rsquared = R2 (estimaciones, prueba$peso),
RMSE = RMSE (estimaciones, prueba$peso),
MAE = MAE (estimaciones, prueba$peso))
resultados$variables <- "Todas"
resultados## Rsquared RMSE MAE variables
## 1 0.8895008 2.674347 2.123038 TodasValidaciones prediciones
# Make predictions
probabilities <- modelo_arboles %>% predict(nltest, type = "response")
predicted.classes <- ifelse(probabilities > 0.1, "1", "0")
# Prediction accuracy
observed.classes <- nltest$peso
mean(predicted.classes == observed.classes, na.rm = TRUE)## [1] 0# Make predictions
probabilities <- modelo_arboles %>% predict(prueba, type = "response")
predicted.classes <- ifelse(probabilities > 0.1, "1", "0")
# Prediction accuracy
observed.classes <- prueba$peso
mean(predicted.classes == observed.classes, na.rm = TRUE)## [1] 0#Metricas con validación cruzada
set.seed(0) #Simulación reproducible
subconjuntos <- trainControl(method = "cv", number = 3) #Repetición de la validación cruzada
modelo_arboles <- train(peso ~ finca + mg + diametro + altura, data=arboles, method="lm", trControl=subconjuntos)
resultados <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados$variables <- "Todos"
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 TodosConclusiones sección modelo inicial
No evidenciamos multicolinealidad entre las variables cuantitativas dado que los porcentajes de corelación son distintos,
El presente modelo contiene multiple variables que arrojan significancia, en especial, de las variables categoricas fincFinca_3 y de las variables cuantitativas altura,
No se cumple 2 de los supuestos de normalidad y no autocorrelación, sin embargo, el modelo presenta muy buen rendimiento tanto en su R2 cuadrado como en su RMSE,
La cantidad de datos evidencia una gran diferencia en el rendimiento cuando aplicamos validación cruzada,
Se explorará la posibilidad de mejorar el modelo realizando una selección de variables, que nos permita tener un modelo más simple y con mayor rendimiento. Dentro de este proceso incluiremos validación cruzada para garantizar que el modelo es el mejor, no omite datos importantes y disminuye el riesgo de sesgos y variabilidad.
Selección de variables
La selección de características las haremos con enfoque en el coeficiente de determinación y el root mean square error (RMSE) despues de un proceso de validación cruzada de 3 iteraciones. Los resultados deberían mejorar los resultados obtenidos en el modelo inicial, en caso contrario, se mantendrían todas las variables.
Altura + Finca
Se propone el presente ejercicio porque son las variables con mayor significacia del modelo.
modelo_arboles <- train(peso ~ finca + altura, data=arboles, method="lm", trControl=subconjuntos)
resultados_to_add <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados_to_add$variables <- "Altura, Finca"
resultados <- rbind(resultados, resultados_to_add)
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 Todos
## 2 0.8087684 3.602604 2.763904 Altura, Finca#No se evidencia una mejora del modelo.
######### Altura + Diametro
#Se propone tener las 2 variables cuantitativas del modelo
modelo_arboles <- train(peso ~ altura + diametro, data=arboles, method="lm", trControl=subconjuntos)
resultados_to_add <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados_to_add$variables <- "Altura, Diametro"
resultados <- rbind(resultados, resultados_to_add)
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 Todos
## 2 0.8087684 3.602604 2.763904 Altura, Finca
## 3 0.8441536 3.628982 2.846408 Altura, Diametro# Se evidencia una mejoría en el modelo.
########### Diametro + Finca
# Como lo observamos en la primera gráfica de esta exploración, el diametro y las fincas formaban un canal muy claro parecido a una linea recta, esto sumado a que el diametro es la variable con mayor correlación y finca es la de mayor significancia podría obtener buen resultado.
modelo_arboles <- train(peso ~ finca + diametro, data=arboles, method="lm", trControl=subconjuntos)
resultados_to_add <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados_to_add$variables <- "Diametro, Finca"
resultados <- rbind(resultados, resultados_to_add)
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 Todos
## 2 0.8087684 3.602604 2.763904 Altura, Finca
## 3 0.8441536 3.628982 2.846408 Altura, Diametro
## 4 0.8523773 3.194301 2.504711 Diametro, Finca# El presente resultado evidencia una aproximación a tener el mismo resultado que todas la variables, estas 2 varibles serán de gran relevancia apartir de ahora.
######## Finca + Diametro + Altura
modelo_arboles <- train(peso ~ finca + diametro + altura, data=arboles, method="lm", trControl=subconjuntos)
resultados_to_add <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados_to_add$variables <- "Finca, Diametro, Altura"
resultados <- rbind(resultados, resultados_to_add)
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 Todos
## 2 0.8087684 3.602604 2.763904 Altura, Finca
## 3 0.8441536 3.628982 2.846408 Altura, Diametro
## 4 0.8523773 3.194301 2.504711 Diametro, Finca
## 5 0.8481022 3.231197 2.457689 Finca, Diametro, Altura# Se disminuye el rendimiento al añadir la altura, finalmente revisaremos el mg.
########### Finca + Diametro + MG
modelo_arboles <- train(peso ~ finca + mg + diametro, data=arboles, method="lm", trControl=subconjuntos)
resultados_to_add <- subset(modelo_arboles$results, select=c(Rsquared, RMSE, MAE))
resultados_to_add$variables <- "Finca, Diametro, MG"
resultados <- rbind(resultados, resultados_to_add)
resultados## Rsquared RMSE MAE variables
## 1 0.8665193 3.127187 2.389579 Todos
## 2 0.8087684 3.602604 2.763904 Altura, Finca
## 3 0.8441536 3.628982 2.846408 Altura, Diametro
## 4 0.8523773 3.194301 2.504711 Diametro, Finca
## 5 0.8481022 3.231197 2.457689 Finca, Diametro, Altura
## 6 0.8455112 3.278228 2.569247 Finca, Diametro, MGDado los resultados obtenidos hasta el momento, optaría por seleccionar la variable finca + diametro, las cuales simplifican el modelo y obtenien un rendimiento muy cercano despues de una validación cruzada de 3 iteracciones al rendimiento de todas las variables, adiciona esto se podría involucrar la altura pero esta teniendo una desmejora pequeña en en R2 y el RMSE. El modelo final sería el siguiente:
modelo_arboles2 <- lm(peso ~ finca + diametro, data=arboles)
summary(modelo_arboles2)##
## Call:
## lm(formula = peso ~ finca + diametro, data = arboles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8658 -2.1882 -0.5434 2.0327 12.7402
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.7033 1.3822 -7.020 4.83e-10 ***
## fincaFINCA_2 1.0459 0.9652 1.084 0.281557
## fincaFINCA_3 3.0739 0.8617 3.567 0.000592 ***
## diametro 4.9758 0.2730 18.225 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.225 on 86 degrees of freedom
## Multiple R-squared: 0.8489, Adjusted R-squared: 0.8437
## F-statistic: 161.1 on 3 and 86 DF, p-value: < 2.2e-16residuales_arboles <- modelo_arboles2$residuals
shapiro.test(residuales_arboles)##
## Shapiro-Wilk normality test
##
## data: residuales_arboles
## W = 0.92047, p-value = 3.785e-05# No cumplimos con el supuesto de normalidad porque nuestro valor p es 8.795e-06, indica que hay significancia para rechaza la hipotesis nula equivalente a normal.
# Varianza constante
lmtest::bptest(modelo_arboles2) ##
## studentized Breusch-Pagan test
##
## data: modelo_arboles2
## BP = 10.961, df = 3, p-value = 0.01194# Cumplimos el supuesto de homocedasticidad porque nuestro valor p es 0.06375, indica que hay noy significancia para rechaza la hipotesis nula equivalente a varianza constante, esto es validado por la gráfica inicial donde no se nota una apertura del canal de los puntos.
# Errores independientes
lmtest::dwtest(modelo_arboles2) ##
## Durbin-Watson test
##
## data: modelo_arboles2
## DW = 1.2304, p-value = 2.333e-05
## alternative hypothesis: true autocorrelation is greater than 0No cumplimos con el supuesto de errores independientes porque nuestro valor p es equivalente a 0.003, indica que hay evidencia para rechaza la hipotesis nula equivalente a no autocorrelación.
Conclusiones
No cumplimos con los supuestos de la regresión lineal y es limitada las inferencias que podemos realizar sobre los coeficientes. Sin embargo, metricas de rendimiento nos indican que los valores estimados son muy aproximados a los reales. Para furturas mejoras, se recomendaría aumentar el tamaño del conjunto de datos, aplicar transformaciones para aproximarse a una varianza constante (que en el primer caso, habia cumplido el susupuesto), cambiar la variable finca a dummies y/o probar otros modelos diferentes a regresión lineal.