Taller Regresión Lineal
Katherin Quiñones % Lizeth Suarez
2022-08-10
Taller de Regresión Lineal
Punto 1.
Se busca realizar la predicción del precio de las acciones de Ecopetrol a partir de la variación del precio del barril de petróleo WTI producido en colombia. Para esto se plantea realizar una regresión lineal.
Se van a considerar los siguietnes datos para su construcción.
datos <- read.csv("C:/Users/lifes/Documents/5. MAESTRIA CIENCIA DE DATOS/2. SEGUNDO SEMESTRE/2. METODOS ESTADISTICOS/1. MODDULO 1/Taller Regresion Lineal/DatosTaller.csv", sep=";", dec = ",")
head(datos)## ï..Fecha PrecioAccionesEcopetrol PrecioPetroleoWTI
## 1 dic-14-15 1090 35.62
## 2 dic-15-15 1170 36.31
## 3 dic-16-15 1160 37.35
## 4 dic-18-15 1230 34.95
## 5 dic-21-15 1155 34.53
## 6 dic-22-15 1165 35.81Inicialmente se analizara de manera descriptiva el comportamiento de ambas variables
m = round(mean(datos[,"PrecioAccionesEcopetrol"]),2)
sd = round(sd(datos[,"PrecioAccionesEcopetrol"]),2)
print(paste("PrecioAccionesEcopetrol", "Media =",m,"y sd =", sd))## [1] "PrecioAccionesEcopetrol Media = 1108.39 y sd = 78.42"Encontramos que en promedio las acciones de Ecopetrol cuestan $1.108 con una desviación estandar de $78.42.
m = round(mean(datos[,"PrecioPetroleoWTI"]),2)
sd = round(sd(datos[,"PrecioPetroleoWTI"]),2)
print(paste("PrecioPetroleoWTI", "Media =",m,"y sd =", sd))## [1] "PrecioPetroleoWTI Media = 35.53 y sd = 2.12"En cuanto al valor medio de un barril de petróleo esta en 35.53 dolares con una desviación estandar de 2.12 dolares.
Con la finalidad de observar la relación de estas dos variables, se grafica un diagrama de dispesión para ver si es posible explicar la variabilidad del precio de las acciones de Ecopetrol mediante la variación del precio del barril de petróleo.
library(ggplot2)
ggplot(datos, aes(x = PrecioPetroleoWTI, y = PrecioAccionesEcopetrol)) +
geom_point() +
geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Se observa, una relación entre esta dos variables, que inicialmente pareciera no ser necesariamente lineal. Sin embargo al analizar el coeficiente de correlación de Pearson este nos da del 70.74% con lo que podriamos considerar que podemos generar la regresión lineal inicialmente.
cor(datos$PrecioAccionesEcopetrol, datos$PrecioPetroleoWTI)## [1] 0.7074373a. Estimación de la Regresión Lineal
Procedemos a estructurar la regresión de la siguiente forma
\[PrecioAccionesEcopetro = \beta_0 + \beta_1 PrecioPetroleoWTI + \epsilon\]
lm = lm(datos$PrecioAccionesEcopetrol~datos$PrecioPetroleoWTI)
coef(lm)## (Intercept) datos$PrecioPetroleoWTI
## 177.76779 26.19213Generando la estimación de la recta de la siguiente forma:
\[PrecioAccionesEcopetro = 177.77+ 26.19* PrecioPetroleoWTI + \epsilon\]
Esta estimación de la recta presenta un \(R^2 = 50.05\), lo que nos dice que el modelo propuesto solo nos explica el 50% de la variabilidad que presenta el precio de las acciones de ecopetrol.
s = summary(lm)
s$r.squared## [1] 0.5004675b. Significancia de Modelo
Al analizar la significacia del modelo, lo hacemos validando que el \(\hat{\beta}_1\) sea significativo, por tanto la hipotesis a evaluar son las siguientes:
\[H_0: \beta_1 = 0\] \[H_1: \beta_1 \neq 0\]
Al realizar la prueba de significancia observamos que el valor p asociado al coeficiente del precio del barril de petróleo es del 1.02%, lo que nos dice que contrastando este resultado con un nivel de significacia del 5%, no se cuenta con evidencia estadistica suficiente para suponer que el \(\beta_1\) estimado no sea significativo.
summary(lm)##
## Call:
## lm(formula = datos$PrecioAccionesEcopetrol ~ datos$PrecioPetroleoWTI)
##
## 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
## datos$PrecioPetroleoWTI 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.001024c. Interpretacion de los betas
En el caso estudiado vemos que los betas estimados se pueden interpretar de la siguiente forma:
\(\beta_0\): En el caso dado donde el precio del barril de petróleo llegase a ser gratis, el valor de las acciones de Ecopetrol costarían $ 177.77.
\(\beta_1\): Este nos quiere decir que por cada dolar adicional que sume el valor del barril de petróleo el precio de las acciones de Ecopetrol incrementaran en $26.192.
d. Analisis de los supuestos del modelo
Linealidad
datos$yest = predict(lm)
ggplot(datos,aes(x = datos$PrecioPetroleoWTI, y = yest)) +
geom_point() +
geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm$residuals),3)))## Warning: Use of `datos$PrecioPetroleoWTI` is discouraged. Use `PrecioPetroleoWTI` instead.
## Use of `datos$PrecioPetroleoWTI` is discouraged. Use `PrecioPetroleoWTI` instead.## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Independencia
Las hipotensis a contrastar son:
\(H_0:\) Los errores no presentan autocorrelación
\(H_1:\) Los errores presentan autocorrelación
library(lmtest)## Warning: package 'lmtest' was built under R version 4.1.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.1.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericdw = dwtest(lm, alternative = "two.sided", iterations = 1000)
plot(lm$residuals, main = paste("dw P-Value = ",round(dw$p.value,3)))
Varianza Constante
Las hipotensis a contrastar son:
\(H_0:\) Los errores tiene homocedasticidad de varianza
\(H_1:\) Los errores no tiene homocedasticidad de varianza
grupos<-numeric()
grupos[datos$yest<1101]<-1
grupos[datos$yest>=1001 & datos$yest<1101]<-2
grupos[datos$yest>=1101]<-3
table(grupos)## grupos
## 1 2 3
## 2 4 12bt = bartlett.test(lm$residuals~grupos)
plot(datos$yest, lm$residuals,
main = paste("dt P-Value = ",round(bt$p.value,3)))
Normalidad
Las hipotesis a contrastar son:
\(H_0:\) Los errores presentan una distribución normal
\(H_1:\) Los errores no presentan una distribución normal
par(mfrow=c(1,2))
qqnorm(lm$residuals)
qqline(lm$residuals)
hist(lm$residuals)
shapiro.test(lm$residuals)##
## Shapiro-Wilk normality test
##
## data: lm$residuals
## W = 0.89259, p-value = 0.04276e. Validez del Modelo
Analizando el modelo propuesto, observamos que usando unicamente el precio del petroleo no podemos llegar a predecir las acciones dado que el modelo propuesto solo explica el 50% de la informacion el precio de las acciones, además el modelo propuesto no cumple con la normalidad, por tanto, se podría explorar una transformación de la variables con la finalidad de poder ver si este supuesto se puede corregir.
Punto 2.
Se busca realizar la predicción del Salario Minimo (SMMLV) a partir de la variación de la inflación en colombia. Para esto se plantea realizar una regresión lineal.
Se van a considerar los siguientes datos para su construcción.
datos <- read.csv("C:/Users/lifes/Documents/5. MAESTRIA CIENCIA DE DATOS/2. SEGUNDO SEMESTRE/2. METODOS ESTADISTICOS/1. MODDULO 1/Taller Regresion Lineal/DatosTallerPunto2.csv", sep=";", dec = ",")
head(datos)## ï..ANIO INFLACION SMLM
## 1 1999 9.23 236460
## 2 2000 8.75 260100
## 3 2001 7.65 286000
## 4 2002 6.99 309000
## 5 2003 6.49 332000
## 6 2004 5.50 358000Inicialmente se analizara de manera descriptiva el comportamiento de ambas variables
m = round(mean(datos[,"INFLACION"]),2)
sd = round(sd(datos[,"INFLACION"]),2)
print(paste("INFLACION", "Media =",m,"y sd =", sd))## [1] "INFLACION Media = 5.35 y sd = 2.32"Encontramos que en promedio la inflación en Colombia es de 5.35% entre 1999 y 2015 con una desviación estandar de 2.32%.
m = round(mean(datos[,"SMLM"]),2)
sd = round(sd(datos[,"SMLM"]),2)
print(paste("SMLM", "Media =",m,"y sd =", sd))## [1] "SMLM Media = 437078.65 y sd = 129182.57"En cuanto al Salario Minimo en Colombia entre 1999 y 2015, el valor promedio es de $ 437.078.7 con una desviación estandar de $ 129.183.
Con la finalidad de observar la relación de estas dos variables, se realiza un diagrama de dispesión con el fin de ver si es posible explicar la variabilidad del salario minimo en funcion de la inflación anual.
Se observa una relación inversa entre la inflación y el salario minimo en Colombia, es decir, entre mayor es la inflación menor ha sido el salario minimo.
library(ggplot2)
ggplot(datos, aes(x = INFLACION, y = SMLM)) +
geom_point() +
geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Al analizar el coeficiente de correlación de Pearson este nos da -70.87% con lo que podriamos considerar que podemos generar la regresión lineal inicialmente.
cor(datos$INFLACION, datos$SMLM)## [1] -0.7086581Estimación de la Regresión Lineal
Procedemos a estructurar la regresión de la siguiente forma
\[SMLM = \beta_0 + \beta_1 Inflación + \epsilon\]
lm = lm(datos$SMLM~datos$INFLACION)
coef(lm)## (Intercept) datos$INFLACION
## 648485.93 -39489.33Generando la estimación de la recta de la siguiente forma:
\[SMLM = 648.485,93 - 39.489,33*Inflacion + \epsilon\]
Significancia del Modelo
Al analizar la significacia del modelo, lo hacemos validando que el \(\hat{\beta}_1\) sea significativo, por tanto la hipotesis a evaluar son las siguientes:
\[H_0: \beta_1 = 0\] \[H_1: \beta_1 \neq 0\]
Al realizar la prueba de significancia observamos que el valor p asociado al coeficiente de la inflación es de 1.45%, lo que nos dice que contrastando este resultado con un nivel de significacia del 5%, no se cuenta con evidencia estadistica suficiente para suponer que el \(\beta_1\) estimado no sea significativo.
summary(lm)##
## Call:
## lm(formula = datos$SMLM ~ datos$INFLACION)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75463 -63456 -42854 17623 263207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 648486 58947 11.00 1.4e-08 ***
## datos$INFLACION -39489 10151 -3.89 0.00145 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 94130 on 15 degrees of freedom
## Multiple R-squared: 0.5022, Adjusted R-squared: 0.469
## F-statistic: 15.13 on 1 and 15 DF, p-value: 0.00145Coeficiente de Correlacion
Esta estimación de la recta presenta un \(R^2 = 50.22\), lo que nos dice que el modelo propuesto solo nos explica el 50% de la variabilidad que presenta el salario minimo.
s = summary(lm)
s$r.squared## [1] 0.5021963### Interpretacion de los betas
En el caso estudiado vemos que los betas estimados se pueden interpretar de la siguiente forma:
\(\beta_0\): En el caso dado donde la inflación llegase a ser 0, el valor del ingreso minimo sería de $ 648.485,93.
\(\beta_1\): Este nos quiere decir que por cada unidad porcentual adicional que incremente la inflación, el salario minimo disminuirá en $39.489,33.
Analisis de los supuestos del modelo
Linealidad
El modelo cumple con el supuesto de linealidad, dado que los estimaciones realizados por este, se encuentran cercanas a los valores reales obervados.
datos$yest = predict(lm)
ggplot(datos,aes(x = datos$INFLACION, y = yest)) +
geom_point() +
geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm$residuals),3)))## Warning: Use of `datos$INFLACION` is discouraged. Use `INFLACION` instead.
## Use of `datos$INFLACION` is discouraged. Use `INFLACION` instead.## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Independencia
Frente a la independecia, pareciera que existiera alguna correalción en el tiempo dado los resultado de los residuales, sin embargo, se podria realizar una prueba formal para validadr esta hipotesis.
plot(lm$residuals)
Varianza Constante
Frente a la homogeneidad de varianza, observamos que si bien tenemos pocos datos parecietan que esto no presentan una tendencia en los datos, por tanto, el modelo cumple con el supuesto de homogeneidad de varainza.
plot(datos$yest, lm$residuals)
Normalidad
Analizando las graficas del comportamiento de los residuales en el QQ plot y en el histograma, podemos observar que el modelo propuesto no cumple con este supuesto. Por tanto, no tenemos certezas sobre la significancia del modelo propuesto, dado que no podemos aplicar prueba de hipotesis sobre los betas estimados.
par(mfrow=c(1,2))
qqnorm(lm$residuals)
qqline(lm$residuals)
hist(lm$residuals)
Convenencia de usar el modelo
Dado los resultado presentados anterior, el modelo propuesto no lo recomendaria para predecir el valor de SMML, dado que este factor depende de muchas otras cosas más, y no unicamente de la inflacion y esto lo podemos evidenciar en los resultados obtendidos donde el \(R^2\) es del alrededor del 50%. Además de la inconsistencia con el comportamiento normal de las variables, dado que se esperaria que a mayor inflación el salario minimo fuera mayor. Quizas se podría tratar de modelar el incremente del salario minimo en lugar del valor neto.
Punto 3
## Warning: package 'kableExtra' was built under R version 4.1.3## Warning: package 'readxl' was built under R version 4.1.3## Warning: package 'Hmisc' was built under R version 4.1.3## Warning: package 'inspectdf' was built under R version 4.1.3## Warning: package 'VIM' was built under R version 4.1.3## Warning: package 'mice' was built under R version 4.1.3## Warning: package 'BBmisc' was built under R version 4.1.3## Warning: package 'car' was built under R version 4.1.3## Warning: package 'carData' was built under R version 4.1.3## Warning: package 'table1' was built under R version 4.1.3La base acontinuación trabajada contiene información de ofertas de vivienda descargadas del portal Fincaraiz, con variables como el precio de vivienda(millones de pesos COP),el area de la vivienda (m2), entre otras variables de interes.
library(readxl)
#Leyendo la base
datos_vivienda <- read_excel("C:/Users/lifes/Documents/5. MAESTRIA CIENCIA DE DATOS/2. SEGUNDO SEMESTRE/2. METODOS ESTADISTICOS/1. MODDULO 1/Taller Regresion Lineal/Datos_Vivienda.xlsx")
attach(datos_vivienda)
names(datos_vivienda)## [1] "Zona" "piso" "Estrato"
## [4] "precio_millon" "Area_contruida" "parqueaderos"
## [7] "Banos" "Habitaciones" "Tipo"
## [10] "Barrio" "cordenada_longitud" "Cordenada_latitud"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 (publicar el informe final en Rpubs presentando código, resultados e interpretaciones).
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 qué?).
Base antes de hacer fitros
head(datos_vivienda)## # A tibble: 6 x 12
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Zona Sur 2 6 880 237 2 5
## 2 Zona Oeste 2 4 1200 800 3 6
## 3 Zona Sur 3 5 250 86 NA 2
## 4 Zona Sur NA 6 1280 346 4 6
## 5 Zona Sur 2 6 1300 600 4 7
## 6 Zona Sur 3 6 513 160 2 4
## # ... with 5 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>A continuación se miraran algunas descriptivas relacionadas a las variables que se debe hacer filtros.
tabla_tipo = tabla_freq(datos_vivienda$Tipo)
tabla_zona = tabla_freq(datos_vivienda$Zona)
kbl(list(tabla_tipo,tabla_zona),caption = "<center><strong>Tipo de Vivienda y Zona</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
tabla_area = descriptivas(datos_vivienda$Area_contruida)
tabla_precio = descriptivas(datos_vivienda$precio_millon)
kbl(list(tabla_area, tabla_precio ),
caption = "<center><strong>Distribución Area y Precio</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
Filtrando Variables
dim(datos_vivienda)[1]## [1] 8322vars <-c('Tipo', 'Zona', 'Area_contruida', 'precio_millon')
cond <-c('Apartamento', 'Zona Norte', 300, 500)
datos_vivienda_sub <- datos_vivienda %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]< cond[[3]],
.data[[vars[[4]]]]< cond[[4]]
)
dim(datos_vivienda_sub)[1]## [1] 315Inicialmente se cuenta con 8322 registros, después de realizar los 4 filtros sugeridos, quedamos con 315 registros, a continuación se presenta un encabezado del dataframe resultante:
head(datos_vivienda_sub)## # A tibble: 6 x 12
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Zona Norte NA 5 340 106 2 2
## 2 Zona Norte 1 3 135 103 1 3
## 3 Zona Norte NA 5 390 102 2 3
## 4 Zona Norte 11 5 350 137 2 3
## 5 Zona Norte NA 5 430 105 NA 3
## 6 Zona Norte NA 3 110 100 NA 1
## # ... with 5 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>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.
Antes de realizar el análiss exploratorio sugerido se plantea hacer algunas recodificaciones necesarias en las variables de interes.
Revisión de datos Faltantes
datos_vivienda_sub$parqueaderos = as.numeric(datos_vivienda_sub$parqueaderos)## Warning: NAs introduced by coercionapply(is.na(datos_vivienda_sub), 2, sum)## Zona piso Estrato precio_millon
## 0 0 0 0
## Area_contruida parqueaderos Banos Habitaciones
## 0 46 0 0
## Tipo Barrio cordenada_longitud Cordenada_latitud
## 0 0 0 0datos_vivienda_sub$parqueaderos[is.na(datos_vivienda_sub$parqueaderos)] <- 0
apply(is.na(datos_vivienda_sub), 2, sum)## Zona piso Estrato precio_millon
## 0 0 0 0
## Area_contruida parqueaderos Banos Habitaciones
## 0 0 0 0
## Tipo Barrio cordenada_longitud Cordenada_latitud
## 0 0 0 0datos_vivienda_sub$parqueaderos = as.numeric(datos_vivienda_sub$parqueaderos)
datos_vivienda_sub$Estrato = as.factor(datos_vivienda_sub$Estrato)
datos_vivienda_sub$parqueaderos_cat=cut(datos_vivienda_sub$parqueaderos, breaks = c(-1,0,4), include.lowest = F, labels = c("sin_parqueadero","con_parqueadero"))
tabla_parque = tabla_freq(datos_vivienda_sub$parqueaderos_cat)
tabla_estrato = tabla_freq(datos_vivienda_sub$Estrato)
kbl(list(tabla_parque, tabla_estrato),
caption = "<center><strong>Distribución Parqueadero y Estrato</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
par(mfrow = c(2,2))
boxplot(precio_millon, horizontal = TRUE, main = "Boxplot Precio Vivienda",col="#A6CEE3" )
hist(precio_millon,main = "Histograma Precio Vivienda" , ylab = "Frecuencia" , xlab = "Millones de pesos",col="#A6CEE3")
boxplot(Area_contruida, horizontal = TRUE, main = "Boxplot Área Construida",col="#A6CEE3" )
hist(Area_contruida,main = "Histograma Área Construida" , ylab = "Frecuencia" , xlab = "m2",col="#A6CEE3")
Análisis Bivariado
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:Hmisc':
##
## subplot## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutfig1 <- plot_ly(data = datos_vivienda_sub, x =~Area_contruida , y = ~precio_millon ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig1 <- fig1 %>% layout(title = 'Precio Vivienda vs Area Construida',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig1## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modefig2 <- plot_ly(datos_vivienda_sub, y = ~precio_millon, color = ~Estrato, type = "box",colors = "Set1")
fig2 <- fig2 %>% layout(title = 'Precio Vivienda vs Estrato',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig2datos_vivienda_sub$parqueaderos_cat = as.factor(datos_vivienda_sub$parqueaderos_cat)
fig3 <- plot_ly(datos_vivienda_sub, y = ~precio_millon, color = ~parqueaderos_cat, type = "box",colors = "Set1")
fig3 <- fig3 %>% layout(title = 'Precio Vivienda vs Parqueaderos',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig3library('fastDummies')
datos_vivienda_sub1 <- dummy_cols(datos_vivienda_sub, select_columns = 'parqueaderos_cat')
datos_vivienda_sub1 <- dummy_cols(datos_vivienda_sub1, select_columns = 'Estrato')
head(datos_vivienda_sub1)## # A tibble: 6 x 19
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Zona Norte NA 5 340 106 2 2
## 2 Zona Norte 1 3 135 103 1 3
## 3 Zona Norte NA 5 390 102 2 3
## 4 Zona Norte 11 5 350 137 2 3
## 5 Zona Norte NA 5 430 105 0 3
## 6 Zona Norte NA 3 110 100 0 1
## # ... with 12 more variables: Habitaciones <dbl>, Tipo <chr>, Barrio <chr>,
## # cordenada_longitud <dbl>, Cordenada_latitud <dbl>, parqueaderos_cat <fct>,
## # parqueaderos_cat_sin_parqueadero <int>,
## # parqueaderos_cat_con_parqueadero <int>, Estrato_3 <int>, Estrato_4 <int>,
## # Estrato_5 <int>, Estrato_6 <int>### 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).
ml=lm(precio_millon~ Area_contruida+Estrato+parqueaderos_cat,data=datos_vivienda_sub)
summary(ml)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos_cat,
## data = datos_vivienda_sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -329.46 -67.87 -14.19 72.27 638.59
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.7088 61.5235 -0.288 0.7737
## Area_contruida 1.6034 0.1958 8.188 7.09e-15 ***
## Estrato4 77.5480 58.8863 1.317 0.1888
## Estrato5 133.2809 56.5569 2.357 0.0191 *
## Estrato6 311.9427 58.6761 5.316 2.03e-07 ***
## parqueaderos_catcon_parqueadero 30.9780 20.4741 1.513 0.1313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 123 on 309 degrees of freedom
## Multiple R-squared: 0.4167, Adjusted R-squared: 0.4073
## F-statistic: 44.15 on 5 and 309 DF, p-value: < 2.2e-16\(H_0 :\beta_1,.., \beta_6=0\)
\(H_a :\beta_1,.., \beta_6\neq 0\)
Se Puede observar que \(|t_0| > t_{\alpha/2 ,n-2}\), entonces se rechaza \(H_0\), y se dice que \(\beta_1\) es diferente de 0.
Observamos que en el modelo inicial reflejan variables que no son “significativas”; no obstante esto puede ser resultante de un problema de multicolinealidad presentado por las variables X; en donde vemos obligados a presentar y esquematizar otro planteamiento que elimine este hallazgo.
Multiple R-squared: 0.4167, R2 bajo!!
El coeficiente de determinación representa la proporción de la variabilidad de Y que es posible explicar a travez de x. En este caso el modelo construido explica el 42% de las variaciones del precio de las viviendas a travez del area construida, el estrato y el número de parqueaderos.
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).
Interpretación Gráfica
par(mfrow=c(2,2))
plot(ml)
Se puede observar que en el gráfico de residuales versus valores
ajustados hay una presencia de valores extremos lo que no hace que del
todo se cumpla que el valor esperado este cerca de 0 en todos los casos.
En cuanto al QQ plot se puede observar que tiene una tendencia lineal
podría pensar en normalidad de no ser por los valores atipicos
presentados en el extremo.
Pruebas estadísticas para validación de supuestos
- \(E(e_i )= 0\)
mean(residuals(ml))## [1] 7.63763e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
shapiro.test(ml$residuals)##
## Shapiro-Wilk normality test
##
## data: ml$residuals
## W = 0.90077, p-value = 1.659e-13Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- Supuesto de Homocedasticidad de Varianza (Breush Pagan)
\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)
\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)
bptest(ml)##
## studentized Breusch-Pagan test
##
## data: ml
## BP = 129.59, df = 5, p-value < 2.2e-16Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
- Supuesto de Autocorrelación de los errores (Durbin-Watson)
\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)
\(H_a : Existe \, correlación\, entre\, los\, errores\)
dwt(ml, alternative = "two.sided")## lag Autocorrelation D-W Statistic p-value
## 1 0.1293537 1.739065 0.02
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), n se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
En este caso se viola el supuesto de normalidad y autocorrelación, por ende se podría pensar en realizar alguna transformación sobre la variable respuesta o inclunsive sobre la variable explicativa.
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?
A continuación buscamos los coeficientes asociados a los betas generados por los modelos.
#(predict(ml,newdata = list(Area_contruida=100,Estrato=1,parqueaderos_cat=1)))
ml$coefficients## (Intercept) Area_contruida
## -17.708836 1.603372
## Estrato4 Estrato5
## 77.547984 133.280900
## Estrato6 parqueaderos_catcon_parqueadero
## 311.942728 30.977981Así pues la ecuación nos queda de la siguiente manera:
\(\hat {precio\_apto}=-17.708836+1.603372 X_{area} + 77.547984 X_{estrato4} + 133.280900 X_{estrato5} + 311.942728 X_{estrato6} +30.977981 X_{con_parqueadero}\)
precio_pred=-17.708836+1.603372*100+ 77.547984*1 + 133.280900*0 + 311.942728*0 +30.977981*1
precio_pred## [1] 251.1543Con base en el modelo propuesto, la oferta del apartamento en 450 millones no parece ser una buena oferta, ya que el precio sugerido, esta casi por la mitad del valor solicitado. (251 millones de pesos).
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.
dim(datos_vivienda)[1]## [1] 8322vars <-c('Tipo', 'Zona', 'Area_contruida', 'precio_millon','parqueaderos')
cond <-c('Apartamento', 'Zona Norte', 100, 400,1)
datos_vivienda_sub_pred <- datos_vivienda %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]> cond[[3]],
.data[[vars[[4]]]]<= cond[[4]],
.data[[vars[[5]]]]>= cond[[5]],
)
ID=1:dim(datos_vivienda_sub_pred)[1]
datos_vivienda_sub_pred=data.frame(ID,datos_vivienda_sub_pred)
dim(datos_vivienda_sub_pred)[1]## [1] 879head(datos_vivienda_sub_pred)## ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1 1 Zona Norte 2 3 135 56 1 1
## 2 2 Zona Norte NA 5 340 106 2 2
## 3 3 Zona Norte 1 3 135 103 1 3
## 4 4 Zona Norte NA 4 175 77 1 2
## 5 5 Zona Norte 2 3 110 57 1 2
## 6 6 Zona Norte 3 3 155 62 1 2
## Habitaciones Tipo Barrio cordenada_longitud
## 1 3 Apartamento torres de comfandi -76.46745
## 2 3 Apartamento la flora -76.48200
## 3 4 Apartamento calimio norte -76.48347
## 4 3 Apartamento urbanizaciv=n pacara -76.48395
## 5 3 Apartamento puente del comercio -76.48458
## 6 3 Apartamento villa del prado -76.48647
## Cordenada_latitud
## 1 3.40763
## 2 3.43500
## 3 3.48626
## 4 3.45104
## 5 3.46987
## 6 3.46549datos_vivienda_sub_pred$parqueaderos = as.numeric(datos_vivienda_sub_pred$parqueaderos)## Warning: NAs introduced by coerciondatos_vivienda_sub_pred$Estrato = as.factor(datos_vivienda_sub_pred$Estrato)
datos_vivienda_sub_pred$parqueaderos[is.na(datos_vivienda_sub_pred$parqueaderos)] <- 0
datos_vivienda_sub_pred$parqueaderos_cat=cut(datos_vivienda_sub_pred$parqueaderos, breaks = c(-1,0,4), include.lowest = F, labels = c("sin_parqueadero","con_parqueadero"))
library('fastDummies')
datos_vivienda_sub_pred1 <- dummy_cols(datos_vivienda_sub_pred, select_columns = 'parqueaderos_cat')
datos_vivienda_sub_pred1 <- dummy_cols(datos_vivienda_sub_pred1, select_columns = 'Estrato')
head(datos_vivienda_sub_pred1)## ID Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1 1 Zona Norte 2 3 135 56 1 1
## 2 2 Zona Norte NA 5 340 106 2 2
## 3 3 Zona Norte 1 3 135 103 1 3
## 4 4 Zona Norte NA 4 175 77 1 2
## 5 5 Zona Norte 2 3 110 57 1 2
## 6 6 Zona Norte 3 3 155 62 1 2
## Habitaciones Tipo Barrio cordenada_longitud
## 1 3 Apartamento torres de comfandi -76.46745
## 2 3 Apartamento la flora -76.48200
## 3 4 Apartamento calimio norte -76.48347
## 4 3 Apartamento urbanizaciv=n pacara -76.48395
## 5 3 Apartamento puente del comercio -76.48458
## 6 3 Apartamento villa del prado -76.48647
## Cordenada_latitud parqueaderos_cat parqueaderos_cat_sin_parqueadero
## 1 3.40763 con_parqueadero 0
## 2 3.43500 con_parqueadero 0
## 3 3.48626 con_parqueadero 0
## 4 3.45104 con_parqueadero 0
## 5 3.46987 con_parqueadero 0
## 6 3.46549 con_parqueadero 0
## parqueaderos_cat_con_parqueadero Estrato_3 Estrato_4 Estrato_5 Estrato_6
## 1 1 1 0 0 0
## 2 1 0 0 1 0
## 3 1 1 0 0 0
## 4 1 0 1 0 0
## 5 1 1 0 0 0
## 6 1 1 0 0 0datos_vivienda_sub_pred1=datos_vivienda_sub_pred1[,c("ID","precio_millon","Area_contruida","parqueaderos_cat_con_parqueadero","Estrato_4","Estrato_5","Estrato_6",'cordenada_longitud','Cordenada_latitud','Barrio')]
head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 1 1 135 56 1 0
## 2 2 340 106 1 0
## 3 3 135 103 1 0
## 4 4 175 77 1 1
## 5 5 110 57 1 0
## 6 6 155 62 1 0
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 1 0 0 -76.46745 3.40763 torres de comfandi
## 2 1 0 -76.48200 3.43500 la flora
## 3 0 0 -76.48347 3.48626 calimio norte
## 4 0 0 -76.48395 3.45104 urbanizaciv=n pacara
## 5 0 0 -76.48458 3.46987 puente del comercio
## 6 0 0 -76.48647 3.46549 villa del pradovars <-c('Estrato_4','parqueaderos_cat_con_parqueadero','Area_contruida')
cond <-c(1,1,100)
datos_vivienda_sub_pred1 <- datos_vivienda_sub_pred1 %>% filter(
.data[[vars[[1]]]]== cond[[1]],
.data[[vars[[2]]]]== cond[[2]],
.data[[vars[[3]]]]>= cond[[3]]
)
print(dim(datos_vivienda_sub_pred1)[1])## [1] 25head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 1 327 380 123 1 1
## 2 515 350 130 1 1
## 3 524 290 108 1 1
## 4 583 185 104 1 1
## 5 602 265 125 1 1
## 6 631 380 126 1 1
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 1 0 0 -76.51437 3.48618 la flora
## 2 0 0 -76.52100 3.49000 la flora
## 3 0 0 -76.52115 3.48930 la flora
## 4 0 0 -76.52300 3.46400 san vicente
## 5 0 0 -76.52353 3.48157 la flora
## 6 0 0 -76.52432 3.48254 la floradatos_vivienda_sub_pred1['prediccion_precio']=round(-17.708836+1.603372*datos_vivienda_sub_pred1$Area_contruida+ 77.547984*datos_vivienda_sub_pred1$Estrato_4 + 133.280900*0 + 311.942728*0 +30.977981*datos_vivienda_sub_pred1$parqueaderos_cat_con_parqueadero)
datos_vivienda_sub_pred1['dif_precio']=datos_vivienda_sub_pred1['prediccion_precio']- datos_vivienda_sub_pred1['precio_millon']
datos_vivienda_sub_pred1=datos_vivienda_sub_pred1[order(datos_vivienda_sub_pred1$dif_precio, decreasing = TRUE), ]
head(datos_vivienda_sub_pred1)## ID precio_millon Area_contruida parqueaderos_cat_con_parqueadero Estrato_4
## 9 678 300 287.00 1 1
## 21 802 280 173.00 1 1
## 4 583 185 104.00 1 1
## 11 685 190 104.12 1 1
## 7 643 270 152.00 1 1
## 13 694 300 163.00 1 1
## Estrato_5 Estrato_6 cordenada_longitud Cordenada_latitud Barrio
## 9 0 0 -76.52673 3.47907 la campiv±a
## 21 0 0 -76.53362 3.46337 santa monica
## 4 0 0 -76.52300 3.46400 san vicente
## 11 0 0 -76.52707 3.46279 san vicente
## 7 0 0 -76.52515 3.46334 versalles
## 13 0 0 -76.52785 3.46701 san vicente
## prediccion_precio dif_precio
## 9 551 251
## 21 368 88
## 4 258 73
## 11 258 68
## 7 335 65
## 13 352 52require(leaflet)## Loading required package: leafletdatos_vivienda_sub_pred1_top5=datos_vivienda_sub_pred1[1:5,]
leaflet() %>% addCircleMarkers(lng = datos_vivienda_sub_pred1_top5$cordenada_longitud,lat = datos_vivienda_sub_pred1_top5$Cordenada_latitud,radius = 0.3,color = "black",label = datos_vivienda_sub_pred1_top5$ID) %>% addTiles()
|
PUnto 4
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.
library(readxl)
#Leyendo la base
datos_arboles <- read_excel("C:/Users/lifes/Documents/5. MAESTRIA CIENCIA DE DATOS/2. SEGUNDO SEMESTRE/2. METODOS ESTADISTICOS/1. MODDULO 1/Taller Regresion Lineal/data arboles.xlsx",col_types = c("text", "text", "numeric","numeric", "numeric"))## Warning: Coercing text to numeric in C2 / R2C3: '13.73'## Warning: Coercing text to numeric in D2 / R2C4: '4.7'## Warning: Coercing text to numeric in C3 / R3C3: '14.58'## Warning: Coercing text to numeric in D3 / R3C4: '5.3'## Warning: Coercing text to numeric in E3 / R3C5: '5.6'## Warning: Coercing text to numeric in C4 / R4C3: '15.88'## Warning: Coercing text to numeric in D4 / R4C4: '4.8'## Warning: Coercing text to numeric in E4 / R4C5: '5.8'## Warning: Coercing text to numeric in C5 / R5C3: '8.99'## Warning: Coercing text to numeric in D5 / R5C4: '3.2'## Warning: Coercing text to numeric in E5 / R5C5: '4.3'## Warning: Coercing text to numeric in C6 / R6C3: '6.99'## Warning: Coercing text to numeric in D6 / R6C4: '2.2'## Warning: Coercing text to numeric in E6 / R6C5: '3.3'## Warning: Coercing text to numeric in C7 / R7C3: '19.34'## Warning: Coercing text to numeric in D7 / R7C4: '6.3'## Warning: Coercing text to numeric in E7 / R7C5: '7.9'## Warning: Coercing text to numeric in C8 / R8C3: '21.44'## Warning: Coercing text to numeric in D8 / R8C4: '6.6'## Warning: Coercing text to numeric in E8 / R8C5: '8.3'## Warning: Coercing text to numeric in C9 / R9C3: '13.81'## Warning: Coercing text to numeric in D9 / R9C4: '5.3'## Warning: Coercing text to numeric in E9 / R9C5: '7.3'## Warning: Coercing text to numeric in C10 / R10C3: '11.88'## Warning: Coercing text to numeric in D10 / R10C4: '4.9'## Warning: Coercing text to numeric in E10 / R10C5: '6.7'## Warning: Coercing text to numeric in C11 / R11C3: '16.62'## Warning: Coercing text to numeric in D11 / R11C4: '5.9'## Warning: Coercing text to numeric in E11 / R11C5: '7.1'## Warning: Coercing text to numeric in C12 / R12C3: '7.03'## Warning: Coercing text to numeric in D12 / R12C4: '2.5'## Warning: Coercing text to numeric in E12 / R12C5: '3.7'## Warning: Coercing text to numeric in C13 / R13C3: '15.83'## Warning: Coercing text to numeric in D13 / R13C4: '4.7'## Warning: Coercing text to numeric in E13 / R13C5: '5.7'## Warning: Coercing text to numeric in C14 / R14C3: '7.47'## Warning: Coercing text to numeric in D14 / R14C4: '2.2'## Warning: Coercing text to numeric in E14 / R14C5: '3.5'## Warning: Coercing text to numeric in C15 / R15C3: '7.87'## Warning: Coercing text to numeric in D15 / R15C4: '3.1'## Warning: Coercing text to numeric in C16 / R16C3: '11.04'## Warning: Coercing text to numeric in D16 / R16C4: '3.7'## Warning: Coercing text to numeric in E16 / R16C5: '4.6'## Warning: Coercing text to numeric in C17 / R17C3: '20.41'## Warning: Coercing text to numeric in D17 / R17C4: '6.5'## Warning: Coercing text to numeric in E17 / R17C5: '7.6'## Warning: Coercing text to numeric in C18 / R18C3: '20.06'## Warning: Coercing text to numeric in D18 / R18C4: '6.4'## Warning: Coercing text to numeric in C19 / R19C3: '18.69'## Warning: Coercing text to numeric in D19 / R19C4: '6.3'## Warning: Coercing text to numeric in E19 / R19C5: '8.1'## Warning: Coercing text to numeric in D20 / R20C4: '4.9'## Warning: Coercing text to numeric in C21 / R21C3: '10.69'## Warning: Coercing text to numeric in D21 / R21C4: '4.6'## Warning: Coercing text to numeric in E21 / R21C5: '6.5'## Warning: Coercing text to numeric in C22 / R22C3: '14.09'## Warning: Coercing text to numeric in D22 / R22C4: '4.2'## Warning: Coercing text to numeric in E22 / R22C5: '5.2'## Warning: Coercing text to numeric in C23 / R23C3: '10.33'## Warning: Coercing text to numeric in D23 / R23C4: '3.6'## Warning: Coercing text to numeric in E23 / R23C5: '4.6'## Warning: Coercing text to numeric in C24 / R24C3: '18.16'## Warning: Coercing text to numeric in D24 / R24C4: '5.6'## Warning: Coercing text to numeric in E24 / R24C5: '6.2'## Warning: Coercing text to numeric in C25 / R25C3: '13.61'## Warning: Coercing text to numeric in D25 / R25C4: '4.6'## Warning: Coercing text to numeric in E25 / R25C5: '5.5'## Warning: Coercing text to numeric in C26 / R26C3: '9.93'## Warning: Coercing text to numeric in D26 / R26C4: '3.5'## Warning: Coercing text to numeric in E26 / R26C5: '4.3'## Warning: Coercing text to numeric in C27 / R27C3: '12.07'## Warning: Coercing text to numeric in E27 / R27C5: '6.4'## Warning: Coercing text to numeric in C28 / R28C3: '9.89'## Warning: Coercing text to numeric in D28 / R28C4: '4.3'## Warning: Coercing text to numeric in E28 / R28C5: '6.3'## Warning: Coercing text to numeric in C29 / R29C3: '6.98'## Warning: Coercing text to numeric in E29 / R29C5: '4.8'## Warning: Coercing text to numeric in C30 / R30C3: '14.44'## Warning: Coercing text to numeric in D30 / R30C4: '5.5'## Warning: Coercing text to numeric in E30 / R30C5: '7.4'## Warning: Coercing text to numeric in C31 / R31C3: '8.73'## Warning: Coercing text to numeric in E31 / R31C5: '5.3'## Warning: Coercing text to numeric in C32 / R32C3: '16.24'## Warning: Coercing text to numeric in D32 / R32C4: '5.8'## Warning: Coercing text to numeric in E32 / R32C5: '8.1'## Warning: Coercing text to numeric in C33 / R33C3: '16.08'## Warning: Coercing text to numeric in D33 / R33C4: '5.9'## Warning: Coercing text to numeric in E33 / R33C5: '7.5'## Warning: Coercing text to numeric in C34 / R34C3: '23.82'## Warning: Coercing text to numeric in D34 / R34C4: '7.1'## Warning: Coercing text to numeric in E34 / R34C5: '9.3'## Warning: Coercing text to numeric in C35 / R35C3: '30.83'## Warning: Coercing text to numeric in D35 / R35C4: '7.9'## Warning: Coercing text to numeric in E35 / R35C5: '10.9'## Warning: Coercing text to numeric in C36 / R36C3: '26.38'## Warning: Coercing text to numeric in D36 / R36C4: '7.1'## Warning: Coercing text to numeric in E36 / R36C5: '9.2'## Warning: Coercing text to numeric in C37 / R37C3: '9.97'## Warning: Coercing text to numeric in D37 / R37C4: '3.7'## Warning: Coercing text to numeric in E37 / R37C5: '4.4'## Warning: Coercing text to numeric in C38 / R38C3: '17.01'## Warning: Coercing text to numeric in D38 / R38C4: '5.4'## Warning: Coercing text to numeric in E38 / R38C5: '6.2'## Warning: Coercing text to numeric in C39 / R39C3: '14.75'## Warning: Coercing text to numeric in D39 / R39C4: '5.1'## Warning: Coercing text to numeric in E39 / R39C5: '5.5'## Warning: Coercing text to numeric in C40 / R40C3: '20.75'## Warning: Coercing text to numeric in D40 / R40C4: '6.2'## Warning: Coercing text to numeric in E40 / R40C5: '6.8'## Warning: Coercing text to numeric in C41 / R41C3: '5.98'## Warning: Coercing text to numeric in D41 / R41C4: '2.5'## Warning: Coercing text to numeric in E41 / R41C5: '3.5'## Warning: Coercing text to numeric in C42 / R42C3: '25.24'## Warning: Coercing text to numeric in D42 / R42C4: '6.6'## Warning: Coercing text to numeric in E42 / R42C5: '8.6'## Warning: Coercing text to numeric in C43 / R43C3: '19.81'## Warning: Coercing text to numeric in D43 / R43C4: '6.9'## Warning: Coercing text to numeric in E43 / R43C5: '8.1'## Warning: Coercing text to numeric in C44 / R44C3: '32.44'## Warning: Coercing text to numeric in D44 / R44C4: '7.9'## Warning: Coercing text to numeric in E44 / R44C5: '9.4'## Warning: Coercing text to numeric in C45 / R45C3: '32.69'## Warning: Coercing text to numeric in D45 / R45C4: '7.8'## Warning: Coercing text to numeric in C46 / R46C3: '16.94'## Warning: Coercing text to numeric in D46 / R46C4: '5.4'## Warning: Coercing text to numeric in C47 / R47C3: '19.02'## Warning: Coercing text to numeric in C48 / R48C3: '14.99'## Warning: Coercing text to numeric in D48 / R48C4: '5.2'## Warning: Coercing text to numeric in E48 / R48C5: '5.8'## Warning: Coercing text to numeric in C49 / R49C3: '22.01'## Warning: Coercing text to numeric in D49 / R49C4: '6.5'## Warning: Coercing text to numeric in C50 / R50C3: '20.24'## Warning: Coercing text to numeric in D50 / R50C4: '6.4'## Warning: Coercing text to numeric in E50 / R50C5: '7.4'## Warning: Coercing text to numeric in C51 / R51C3: '20.06'## Warning: Coercing text to numeric in E51 / R51C5: '6.6'## Warning: Coercing text to numeric in C52 / R52C3: '30.51'## Warning: Coercing text to numeric in D52 / R52C4: '7.2'## Warning: Coercing text to numeric in E52 / R52C5: '9.4'## Warning: Coercing text to numeric in C53 / R53C3: '33.42'## Warning: Coercing text to numeric in D53 / R53C4: '7.7'## Warning: Coercing text to numeric in E53 / R53C5: '10.1'## Warning: Coercing text to numeric in C54 / R54C3: '45.41'## Warning: Coercing text to numeric in D54 / R54C4: '8.7'## Warning: Coercing text to numeric in E54 / R54C5: '9.6'## Warning: Coercing text to numeric in C55 / R55C3: '37.74'## Warning: Coercing text to numeric in D55 / R55C4: '8.2'## Warning: Coercing text to numeric in E55 / R55C5: '10.5'## Warning: Coercing text to numeric in C56 / R56C3: '47.87'## Warning: Coercing text to numeric in D56 / R56C4: '8.8'## Warning: Coercing text to numeric in E56 / R56C5: '11.3'## Warning: Coercing text to numeric in C57 / R57C3: '12.95'## Warning: Coercing text to numeric in D57 / R57C4: '4.8'## Warning: Coercing text to numeric in E57 / R57C5: '5.5'## Warning: Coercing text to numeric in C58 / R58C3: '17.71'## Warning: Coercing text to numeric in D58 / R58C4: '5.8'## Warning: Coercing text to numeric in C59 / R59C3: '22.75'## Warning: Coercing text to numeric in D59 / R59C4: '6.6'## Warning: Coercing text to numeric in E59 / R59C5: '6.1'## Warning: Coercing text to numeric in C60 / R60C3: '23.02'## Warning: Coercing text to numeric in D60 / R60C4: '6.8'## Warning: Coercing text to numeric in E60 / R60C5: '7.5'## Warning: Coercing text to numeric in C61 / R61C3: '14.03'## Warning: Coercing text to numeric in E61 / R61C5: '5.7'## Warning: Coercing text to numeric in C62 / R62C3: '13.98'## Warning: Coercing text to numeric in D62 / R62C4: '4.7'## Warning: Coercing text to numeric in E62 / R62C5: '5.8'## Warning: Coercing text to numeric in C63 / R63C3: '24.47'## Warning: Coercing text to numeric in E63 / R63C5: '7.2'## Warning: Coercing text to numeric in C64 / R64C3: '22.82'## Warning: Coercing text to numeric in D64 / R64C4: '6.7'## Warning: Coercing text to numeric in E64 / R64C5: '7.8'## Warning: Coercing text to numeric in C65 / R65C3: '21.8'## Warning: Coercing text to numeric in D65 / R65C4: '6.8'## Warning: Coercing text to numeric in E65 / R65C5: '7.3'## Warning: Coercing text to numeric in C66 / R66C3: '18.12'## Warning: Coercing text to numeric in D66 / R66C4: '5.9'## Warning: Coercing text to numeric in E66 / R66C5: '7.7'## Warning: Coercing text to numeric in C67 / R67C3: '22.85'## Warning: Coercing text to numeric in E67 / R67C5: '7.3'## Warning: Coercing text to numeric in C68 / R68C3: '23.14'## Warning: Coercing text to numeric in E68 / R68C5: '7.5'## Warning: Coercing text to numeric in C69 / R69C3: '27.33'## Warning: Coercing text to numeric in D69 / R69C4: '6.4'## Warning: Coercing text to numeric in C70 / R70C3: '29.16'## Warning: Coercing text to numeric in E70 / R70C5: '8.3'## Warning: Coercing text to numeric in C71 / R71C3: '17.26'## Warning: Coercing text to numeric in D71 / R71C4: '5.2'## Warning: Coercing text to numeric in E71 / R71C5: '6.5'## Warning: Coercing text to numeric in C72 / R72C3: '34.38'## Warning: Coercing text to numeric in E72 / R72C5: '8.8'## Warning: Coercing text to numeric in C73 / R73C3: '27.87'## Warning: Coercing text to numeric in D73 / R73C4: '6.7'## Warning: Coercing text to numeric in E73 / R73C5: '8.3'## Warning: Coercing text to numeric in C74 / R74C3: '27.45'## Warning: Coercing text to numeric in D74 / R74C4: '6.8'## Warning: Coercing text to numeric in E74 / R74C5: '7.2'## Warning: Coercing text to numeric in C75 / R75C3: '31.16'## Warning: Coercing text to numeric in D75 / R75C4: '6.7'## Warning: Coercing text to numeric in E75 / R75C5: '8.5'## Warning: Coercing text to numeric in C76 / R76C3: '23.35'## Warning: Coercing text to numeric in D76 / R76C4: '5.9'## Warning: Coercing text to numeric in E76 / R76C5: '7.5'## Warning: Coercing text to numeric in C77 / R77C3: '14.52'## Warning: Coercing text to numeric in D77 / R77C4: '4.4'## Warning: Coercing text to numeric in E77 / R77C5: '4.9'## Warning: Coercing text to numeric in C78 / R78C3: '13.01'## Warning: Coercing text to numeric in D78 / R78C4: '4.2'## Warning: Coercing text to numeric in E78 / R78C5: '4.8'## Warning: Coercing text to numeric in C79 / R79C3: '16.19'## Warning: Coercing text to numeric in D79 / R79C4: '4.7'## Warning: Coercing text to numeric in E79 / R79C5: '5.6'## Warning: Coercing text to numeric in C80 / R80C3: '18.67'## Warning: Coercing text to numeric in D80 / R80C4: '5.3'## Warning: Coercing text to numeric in C81 / R81C3: '13.23'## Warning: Coercing text to numeric in D81 / R81C4: '4.4'## Warning: Coercing text to numeric in E81 / R81C5: '4.6'## Warning: Coercing text to numeric in C82 / R82C3: '14.03'## Warning: Coercing text to numeric in D82 / R82C4: '4.2'## Warning: Coercing text to numeric in E82 / R82C5: '4.7'## Warning: Coercing text to numeric in C83 / R83C3: '13.3'## Warning: Coercing text to numeric in D83 / R83C4: '3.8'## Warning: Coercing text to numeric in E83 / R83C5: '4.3'## Warning: Coercing text to numeric in C84 / R84C3: '17.96'## Warning: Coercing text to numeric in D84 / R84C4: '5.3'## Warning: Coercing text to numeric in E84 / R84C5: '5.8'## Warning: Coercing text to numeric in C85 / R85C3: '18.67'## Warning: Coercing text to numeric in D85 / R85C4: '5.2'## Warning: Coercing text to numeric in C86 / R86C3: '13.34'## Warning: Coercing text to numeric in D86 / R86C4: '3.5'## Warning: Coercing text to numeric in C87 / R87C3: '15.4'## Warning: Coercing text to numeric in D87 / R87C4: '4.1'## Warning: Coercing text to numeric in E87 / R87C5: '4.6'## Warning: Coercing text to numeric in C88 / R88C3: '14.91'## Warning: Coercing text to numeric in D88 / R88C4: '3.9'## Warning: Coercing text to numeric in E88 / R88C5: '4.8'## Warning: Coercing text to numeric in C89 / R89C3: '17.73'## Warning: Coercing text to numeric in E89 / R89C5: '5.2'## Warning: Coercing text to numeric in C90 / R90C3: '21.12'## Warning: Coercing text to numeric in D90 / R90C4: '5.4'## Warning: Coercing text to numeric in C91 / R91C3: '18.49'## Warning: Coercing text to numeric in D91 / R91C4: '4.5'## Warning: Coercing text to numeric in E91 / R91C5: '5.1'attach(datos_arboles)
print(names(datos_arboles))## [1] "finca" "mg" "peso" "diametro" "altura"head(datos_arboles)## # A tibble: 6 x 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.9Se puede observar que en el archivo dispuesto, se encuentran variables como la finca, el tipo de genotipo, estas como variables cualitativas y como variables cuantitativas las variables peso, diametro, y altura de los arboles.. A continuaciòn se hara la construcciòn de dos modelos de regresiòn simple usando las variables continua y concluiremos cual podrìa ser el mejor modelo.
Lo primero que se revisarà es si la base de datos contiene datos faltantes y de ser asì, se buscara hacer un tratamiento de ellos, previo a la construcciòn del ejercicio propuesto.
apply(is.na(datos_arboles), 2, sum)## finca mg peso diametro altura
## 0 0 0 0 0tabla_finca = tabla_freq(datos_arboles$finca)
tabla_gen = tabla_freq(datos_arboles$mg)
kbl(list(tabla_finca,tabla_gen),caption = "<center><strong>Tipo de finca y Genotipo árbol</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
datos_arboles$peso=as.numeric(datos_arboles$peso)
datos_arboles$diametro=as.numeric(datos_arboles$diametro)
datos_arboles$altura=as.numeric(datos_arboles$altura)
tabla_peso = descriptivas(datos_arboles$peso)
tabla_diametro = descriptivas(datos_arboles$diametro)
tabla_altura = descriptivas(datos_arboles$altura)
kbl(list(tabla_peso, tabla_diametro,tabla_altura ),
caption = "<center><strong>Distribución Peso, Diametro, Altura de Àrboles</strong></center>") %>%
kable_paper(bootstrap_options = "striped")
|
|
|
Análisis Bivariado
library(plotly)
fig1 <- plot_ly(data = datos_arboles, x =~altura , y = ~peso ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig1 <- fig1 %>% layout(title = 'Peso vs Altura de Àrbol',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig1## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modelibrary(plotly)
fig2 <- plot_ly(data = datos_arboles, x =~diametro , y = ~peso ,
marker = list(size = 10,
color = 'lightblue',
line = list(color = 'blue',
width = 2)))
fig2 <- fig2 %>% layout(title = 'Peso vs Diametro de Arbol',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
fig2## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modeConstrucción Modelos:
Modelo 1: Altura vs Peso
m1=lm(peso ~altura,data=datos_arboles)
summary(m1)##
## Call:
## lm(formula = peso ~ altura, data = datos_arboles)
##
## 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 ***
## altura 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-16cor1=cor(datos_arboles$peso,datos_arboles$altura)
print('Correlación Peso vs Altura:',str(cor1))## num 0.858
## [1] "Correlación Peso vs Altura:"rcuad_m1=cor1*cor1
print('Coeficiente de Determinacion Peso vs Diametro:',str(rcuad_m1))## num 0.737
## [1] "Coeficiente de Determinacion Peso vs Diametro:"Modelo 2: Diametro vs Peso
m2=lm(peso ~diametro,data=datos_arboles)
summary(m2)##
## Call:
## lm(formula = peso ~ diametro, data = datos_arboles)
##
## 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 ***
## diametro 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-16cor2=cor(datos_arboles$peso,datos_arboles$diametro)
rcuad_m2=cor2*cor2
print('Correlación Peso vs Diametro:',str(cor2))## num 0.908
## [1] "Correlación Peso vs Diametro:"print('Coeficiente de Determinacion Peso vs Diametro:',str(rcuad_m2))## num 0.825
## [1] "Coeficiente de Determinacion Peso vs Diametro:"Modelo 3: Diametro + Altura vs Peso
m3=lm(peso ~diametro + altura ,data=datos_arboles)
summary(m3)##
## Call:
## lm(formula = peso ~ diametro + altura, data = datos_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 ***
## diametro 4.7395 0.7128 6.649 2.49e-09 ***
## altura 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-16Validación de Supuestos Modelo 1:
- \(E(e_i )= 0\)
error_medio_m1=mean(residuals(m1))
error_medio_m1## [1] -1.28331e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m1=shapiro.test(m1$residuals)
norm_m1##
## Shapiro-Wilk normality test
##
## data: m1$residuals
## W = 0.95439, p-value = 0.003157Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- Supuesto de Homocedasticidad de Varianza (Breush Pagan)
\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)
\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)
homocedasticidad_m1=bptest(m1)
homocedasticidad_m1##
## studentized Breusch-Pagan test
##
## data: m1
## BP = 12.716, df = 1, p-value = 0.0003625Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
- Supuesto de Autocorrelación de los errores (Durbin-Watson)
\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)
\(H_a : Existe \, correlación\, entre\, los\, errores\)
autocor_m1=dwt(m1, alternative = "two.sided")
autocor_m1## lag Autocorrelation D-W Statistic p-value
## 1 0.537971 0.9021616 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Validación de Supuestos Modelo 2:
- \(E(e_i )= 0\)
error_medio_m2=mean(residuals(m2))
error_medio_m2## [1] 2.553898e-17El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m2=shapiro.test(m2$residuals)
norm_m2##
## Shapiro-Wilk normality test
##
## data: m2$residuals
## W = 0.95356, p-value = 0.002793Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- Supuesto de Homocedasticidad de Varianza (Breush Pagan)
\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)
\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)
homocedasticidad_m2=bptest(m2)
homocedasticidad_m2##
## studentized Breusch-Pagan test
##
## data: m2
## BP = 11.76, df = 1, p-value = 0.0006052Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
- Supuesto de Autocorrelación de los errores (Durbin-Watson)
\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)
\(H_a : Existe \, correlación\, entre\, los\, errores\)
autocor_m2=dwt(m2, alternative = "two.sided")
autocor_m2## lag Autocorrelation D-W Statistic p-value
## 1 0.4533746 1.071861 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Validación de Supuestos Modelo 3:
- \(E(e_i )= 0\)
error_medio_m3=mean(residuals(m3))
error_medio_m3## [1] 1.032498e-16El valor esperado (promedio) de los residuales es igual a 0.
- Supuesto de Normalidad
\(H_0 : X \sim N(\mu,\sigma^2)\)
\(H_a : X \nsim N(\mu,\sigma^2)\)
norm_m3=shapiro.test(m3$residuals)
norm_m3##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.95745, p-value = 0.004966Com \(P-value\) es menor a a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces podría pensar que los errores no siguen una distribución normal.
- Supuesto de Homocedasticidad de Varianza (Breush Pagan)
\(H_0 : Los \, residuales \, se \, distribuyen \, con \, la \, misma \, varianza\)
\(H_a : Los \, residuales \, no \ se \, distribuyen \, con \, la \, misma \, varianza\)
homocedasticidad_m3=bptest(m3)
homocedasticidad_m3##
## studentized Breusch-Pagan test
##
## data: m3
## BP = 14.32, df = 2, p-value = 0.0007772Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza \(H_0\), entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
- Supuesto de Autocorrelación de los errores (Durbin-Watson)
\(H_0 : No \, existe \, correlación \, entre \, los \, errores\)
\(H_a : Existe \, correlación\, entre\, los\, errores\)
autocor_m3=dwt(m3, alternative = "two.sided")
autocor_m3## lag Autocorrelation D-W Statistic p-value
## 1 0.4648495 1.048132 0
## Alternative hypothesis: rho != 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza \(H_0\), entonces se podría pensar que los errores podrian estar autocorrelacionados.
Revisión Criterio Akaike y Bayesiano
akaike_m1=AIC(m1)
bayesian_m1=BIC(m1)
print('Criterio Akaike m1',str(akaike_m1))## num 518
## [1] "Criterio Akaike m1"print('Criterio Bayesian m1',str(bayesian_m1))## num 526
## [1] "Criterio Bayesian m1"akaike_m2=AIC(m2)
bayesian_m2=BIC(m2)
print('Criterio Akaike m2',str(akaike_m2))## num 482
## [1] "Criterio Akaike m2"print('Criterio Bayesian m2',str(bayesian_m2))## num 489
## [1] "Criterio Bayesian m2"akaike_m3=AIC(m3)
bayesian_m3=BIC(m3)
print('Criterio Akaike m3',str(akaike_m3))## num 483
## [1] "Criterio Akaike m3"print('Criterio Bayesian m3',str(bayesian_m3))## num 493
## [1] "Criterio Bayesian m3"Generacion de valores aleatorios para prediccion
y=runif(90,min(datos_arboles$peso),max(datos_arboles$peso))
x_altura=runif(90,min(datos_arboles$altura),max(datos_arboles$altura))
x_diametro=runif(90,min(datos_arboles$diametro),max(datos_arboles$diametro))
datos_arboles2=data.frame(peso=y,altura=x_altura)
head(datos_arboles2)## peso altura
## 1 14.285413 8.677585
## 2 8.271822 6.213213
## 3 7.069356 9.279065
## 4 28.051459 10.126319
## 5 36.528710 7.407937
## 6 23.820860 10.041718Haciendo Prediccion Modelo 1
prediccion_m1=predict(m1,newdata = datos_arboles2 )
resultados_m1=data.frame(y_real=y,y_predict_m1=prediccion_m1)
head(resultados_m1)## y_real y_predict_m1
## 1 14.285413 26.71506
## 2 8.271822 17.12729
## 3 7.069356 29.05516
## 4 28.051459 32.35144
## 5 36.528710 21.77543
## 6 23.820860 32.02230library(Metrics)
attach(resultados_m1)
SCE_m1=sum((y_predict_m1-y)^2)
print('SCE Modelo 1',str(SCE_m1))## num 16360
## [1] "SCE Modelo 1"MAE_m1=mae(y_real,y_predict_m1)
print('MAE Modelo 1',str(MAE_m1))## num 10.5
## [1] "MAE Modelo 1"RMSE_m1=rmse(y_real,y_predict_m1)
print('RMSE Modelo 1',str(RMSE_m1))## num 13.5
## [1] "RMSE Modelo 1"Haciendo Prediccion Modelo 2
datos_arboles3=data.frame(peso=y,diametro=x_diametro)
head(datos_arboles3)## peso diametro
## 1 14.285413 5.052913
## 2 8.271822 3.393959
## 3 7.069356 6.460802
## 4 28.051459 7.977607
## 5 36.528710 7.298708
## 6 23.820860 8.039213prediccion_m2=predict(m2,newdata = datos_arboles3 )
resultados_m2=data.frame(y_real=y,y_predict_m2=prediccion_m2)
head(resultados_m2)## y_real y_predict_m2
## 1 14.285413 16.762621
## 2 8.271822 8.297681
## 3 7.069356 23.946487
## 4 28.051459 31.686099
## 5 36.528710 28.221966
## 6 23.820860 32.000449attach(resultados_m2)## The following object is masked from resultados_m1:
##
## y_realSCE_m2=sum((y_predict_m2-y)^2)
print('SCE Modelo 2',str(SCE_m2))## num 24144
## [1] "SCE Modelo 2"MAE_m2=mae(resultados_m2$y_real,resultados_m2$y_predict_m2)
print('MAE Modelo 2',str(MAE_m2))## num 13.6
## [1] "MAE Modelo 2"RMSE_m2=rmse(resultados_m2$y_real,resultados_m2$y_predict_m2)
print('RMSE Modelo 2',str(RMSE_m2))## num 16.4
## [1] "RMSE Modelo 2"Haciendo Prediccion Modelo 3
datos_arboles4=data.frame(peso=y,diametro=x_diametro,altura=x_altura)
head(datos_arboles4)## peso diametro altura
## 1 14.285413 5.052913 8.677585
## 2 8.271822 3.393959 6.213213
## 3 7.069356 6.460802 9.279065
## 4 28.051459 7.977607 10.126319
## 5 36.528710 7.298708 7.407937
## 6 23.820860 8.039213 10.041718prediccion_m3=predict(m3,newdata = datos_arboles4)
resultados_m3=data.frame(y_real=datos_arboles4$peso,y_predict_m3=prediccion_m3)
head(resultados_m3)## y_real y_predict_m3
## 1 14.285413 17.545029
## 2 8.271822 8.910744
## 3 7.069356 24.406022
## 4 28.051459 31.860179
## 5 36.528710 27.791276
## 6 23.820860 32.125665attach(resultados_m3)## The following object is masked from resultados_m2:
##
## y_real## The following object is masked from resultados_m1:
##
## y_realMAE_m3=mae(resultados_m3$y_real,resultados_m3$y_predict_m3)
print('MAE Modelo 3',str(MAE_m3))## num 13.2
## [1] "MAE Modelo 3"RMSE_m3=rmse(resultados_m3$y_real,resultados_m2$y_predict_m3)
print('RMSE Modelo 3',str(RMSE_m3))## num NaN
## [1] "RMSE Modelo 3"Validación Cruzada Modelo 1
library(rlang)##
## Attaching package: 'rlang'## The following object is masked from 'package:Metrics':
##
## ll## The following object is masked from 'package:tinytex':
##
## check_installed## The following object is masked from 'package:BBmisc':
##
## names2library(vctrs)##
## Attaching package: 'vctrs'## The following object is masked from 'package:dplyr':
##
## data_framelibrary(caret)##
## Attaching package: 'caret'## The following objects are masked from 'package:Metrics':
##
## precision, recall## The following object is masked from 'package:survival':
##
## clustercontrol <- trainControl(method = "cv", number = 5)
#fit a regression model and use k-fold CV to evaluate performance
model_cv_1 <- train(peso ~ altura, data = datos_arboles, method = "lm", trControl = control)
#view summary of k-fold CV
print(model_cv_1)## Linear Regression
##
## 90 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 73, 72, 72, 72, 71
## Resampling results:
##
## RMSE Rsquared MAE
## 4.26194 0.7576605 3.157514
##
## Tuning parameter 'intercept' was held constant at a value of TRUEValidación Cruzada Modelo 2
control <- trainControl(method = "cv", number = 5)
#fit a regression model and use k-fold CV to evaluate performance
model_cv_2 <- train(peso ~ diametro, data = datos_arboles, method = "lm", trControl = control)
#view summary of k-fold CV
print(model_cv_2)## Linear Regression
##
## 90 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 72, 72, 71, 73, 72
## Resampling results:
##
## RMSE Rsquared MAE
## 3.529547 0.8516586 2.881968
##
## Tuning parameter 'intercept' was held constant at a value of TRUEValidación Cruzada Modelo 3
control <- trainControl(method = "cv", number = 5)
#fit a regression model and use k-fold CV to evaluate performance
model_cv_3 <- train(peso ~ diametro + altura, data = datos_arboles, method = "lm", trControl = control)
#view summary of k-fold CV
print(model_cv_3)## Linear Regression
##
## 90 samples
## 2 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 72, 73, 71, 74, 70
## Resampling results:
##
## RMSE Rsquared MAE
## 3.50206 0.8367457 2.891037
##
## Tuning parameter 'intercept' was held constant at a value of TRUEConsolidación de Resultados entre modelos propuestos:
Algunas definiciones de métricas utilizadas se presentan a continuación:
RMSE: La raíz del error cuadrático medio. Esto mide la diferencia promedio entre las predicciones hechas por el modelo y las observaciones reales. Cuanto menor sea el RMSE, más fielmente podrá un modelo predecir las observaciones reales.
Rsquared: Esta es una medida de la correlación entre las predicciones hechas por el modelo y las observaciones reales. Cuanto mayor sea el R-cuadrado, más fielmente podrá un modelo predecir las observaciones reales.
MAE: El error absoluto medio. Esta es la diferencia absoluta promedio entre las predicciones hechas por el modelo y las observaciones reales. Cuanto más bajo es el MAE, más cerca puede un modelo predecir las observaciones reales.
| num_modelo | modelo | correlaciones | coef_det | betas | mae | rmse | mae_cv | rmse_cv | normalidad | heterocedasticidad | autocorrelacion |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Modelo 1 | Altura vs Peso | 0.8582009 | 0.7365088 | Betas Significativos | 10.51297 | 13.48243 | 3.133 | 4.101 | No Cumple | Si Cumple | No Cumple |
| Modelo 2 | Diametro vs Peso | 0.9081230 | 0.8246874 | Betas Significativos | 13.63084 | 16.37891 | 2.769 | 3.377 | No Cumple | Si Cumple | No Cumple |
| Modelo 3 | Diametro + Altura vs Peso | NaN | 0.8253000 | B_2 No Significativo | 13.16665 | NaN | 2.854 | 3.531 | No Cumple | Si Cumple | No Cumple |