Informe de Aplicación de Regresión Lineal
Carlos Andrés Gómez Vasco
2022-08-15
Informe 1 – Modelo de Regresión Lineal Simple y Múltiple
1. Predicción de los precios de las acciones. 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
Cargar Datos
library(readr)
petroleoWTI <- read_delim("D:/MAESTRIA CIENCIA DE DATOS/Segundo Semestre/Modelos Estadisticos/EntregaUnidad1/BasesdeDatos/ptroleoWTI.csv",
delim = ";", escape_double = FALSE, col_types = cols(Precio_Acciones = col_number(),
Precio_WTI= col_number()), trim_ws = TRUE)
petroleoWTI## # A tibble: 18 × 3
## `Fecha ` Precio_Acciones Precio_WTI
## <chr> <dbl> <dbl>
## 1 dic 14-15 1090 35.6
## 2 dic 15-15 1170 36.3
## 3 dic 16-15 1160 37.3
## 4 dic 18-15 1230 34.9
## 5 dic 21-15 1155 34.5
## 6 dic 22-15 1165 35.8
## 7 dic 23-15 1205 36.1
## 8 dic 24-15 1170 37.5
## 9 dic 28-15 1150 37.8
## 10 dic 29-15 1130 36.8
## 11 dic 30-15 1110 37.9
## 12 ene 04-16 1105 37.0
## 13 ene 05-16 1085 36.8
## 14 ene 06-16 1060 36.0
## 15 ene 07-16 1035 34.0
## 16 ene 08-16 1015 33.3
## 17 ene 12-16 955 31.4
## 18 ene 13-16 961 30.4La base de datos muestra el comportamiento en el tiempo de: el precio en pesos de las acciones de Ecopetrol y el respectivo precio del barril en dolares de referencia WTI
Analisis Exploratorio
Al realizar una exploración de los datos se puede observar que el valor promedio de la acción esta en 1108 pesos, mientras el valor promedio del barril es de 35.53 dolares. La distribución de los datos en ambos casos no son muy parecidas dadas las respectivas desviaciones estandar de 78.42 y 2.118183, sin embargi se pueden ilustrar variaciones de las medidas de tendencia central entre un 6% y 7% en ambos casos. Ello podria ser un ligero indicio de la dependencia que podrian tener las variables de Precio_Acciones y Precio_WTI.
attach(petroleoWTI)
summary(petroleoWTI)## Fecha Precio_Acciones Precio_WTI
## Length:18 Min. : 955 Min. :30.44
## Class :character 1st Qu.:1066 1st Qu.:34.63
## Mode :character Median :1120 Median :36.05
## Mean :1108 Mean :35.53
## 3rd Qu.:1164 3rd Qu.:36.98
## Max. :1230 Max. :37.87hist(Precio_Acciones,col="blue")
hist(Precio_WTI,col="blue")
sd(Precio_Acciones)## [1] 78.42354sd(Precio_WTI)## [1] 2.118183require(table1)
table1(~Precio_Acciones+Precio_WTI,data=petroleoWTI)| Overall (N=18) |
|
|---|---|
| Precio_Acciones | |
| Mean (SD) | 1110 (78.4) |
| Median [Min, Max] | 1120 [955, 1230] |
| Precio_WTI | |
| Mean (SD) | 35.5 (2.12) |
| Median [Min, Max] | 36.1 [30.4, 37.9] |
- 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 \[R^{2}\]
Para establecer el modelo primero realizaremos un grafico que nos muestre la tendencia de los datos.
library(ggplot2)
ggplot(petroleoWTI, aes(Precio_WTI, Precio_Acciones)) + geom_point()+scale_x_continuous("Precio Barril Referencia WTI (en dolares)") + scale_y_continuous("Precio Acción Ecopetrol (en pesos)") + geom_smooth()
En el grafico se puede observar que la linea suavizada no muestra una
tendencia de caracter lineal entre los datos. Sin embargo podemos
calcular el coeficiente de correlación entre las dos variables.
cor(petroleoWTI$Precio_WTI,petroleoWTI$Precio_Acciones)## [1] 0.7074373La correlación que se observa no es muy fuerte entre las dos variables, pero tampoco muestra una correlación despreciable. Para mejorar el analisis es preciso establecer el modelo de estimación lineal.
Estimación del Modelo Lineal Simple
En esta estimación se encuentra que la ecuación del moelo ajustado esta dada por: \[PrecioAcciónEcopetrol=177.768+26.192*PrecioBarrilWTI\]
Donde \(R^{2}=0.4692\), es decir; el modelo consigue explicar la correlación de los datos en casi la mitad de los mismos.
modelo1=lm(Precio_Acciones~Precio_WTI,data=petroleoWTI)
summary(modelo1)##
## Call:
## lm(formula = Precio_Acciones ~ Precio_WTI, data = petroleoWTI)
##
## 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
## Precio_WTI 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.001024Una grafica del ajuste en los datos se puede observa como sigue, con un intervalo de confianza del 95% establecido por la sección sombreada en gris:
library(ggplot2)
ggplot(petroleoWTI, aes(Precio_WTI, Precio_Acciones)) + geom_point()+scale_x_continuous("Precio Barril Referencia WTI (en dolares)") + scale_y_continuous("Precio Acción Ecopetrol (en pesos)") + geom_smooth(method = "lm")
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 \(\alpha=0.05\).
De acuerdo con los resultados obtenidos, una hipotesis acorde con el estudio podria estar relacionda con el hecho esencial que si sube el precio del barril de referencia WTI, estaria aumentando de forma lineal el valor de la acción de Ecopetrol. Sin embargo, es evidente que la situación no estan simple y que el modelo lineal no responde a mas de la mitad de los datos para que esta hipotesis sea valida. Aunque de acuerdo al creiterio de \(\alpha=0.05\) y revisando el valor p de significacia que obtenemos en el modelo \(0.001024\) podriamos validar la hipotesis establecida.
- Interprete los coeficientes del modelo propuesto en “a)”
En el modelo conseguido los coeficientes definirian qu para un valor cero del precio del barril WTI la acción de ecopetrol seria equivalente a \(\$177.768\) pesos y adicionalmente por cada dolar que aumente el precio del barril, se consigue aumentar el preciode la acción de ecopetrol en \(\$26.192\) pesos.
- Haga un análisis de los residuos. ¿Qué supuesto no se cumple?
Para validar los supuestos, se realiza el analisis de rersiduos y se encuentra que no cumple con con la media cero aunque si es visible que no hay ninguna tendencia entre los residuales y los ajustados.
Por otra parte existe un ajuste aprobable en la recta de normalidad qqplot. Y al realizar las pruebas de Shapiro para cada variable, el p valor (0.3518) es mayor a \(\alpha=0.05\), no se rechaza que la variable Precio_Acciones presenta un comportamiento normal o paramétrico. Caso contrario ocurre con la variable Precio_WTI con el el p valor 0.03224, que no se comporta de manera normal, confirmando la observación detallada en el analisis estadistico descriptivo de la parte inicial.
par(mfrow=c(2,2))
plot(modelo1)
shapiro.test(petroleoWTI$Precio_Acciones)##
## Shapiro-Wilk normality test
##
## data: petroleoWTI$Precio_Acciones
## W = 0.94498, p-value = 0.3518shapiro.test(petroleoWTI$Precio_WTI)##
## Shapiro-Wilk normality test
##
## data: petroleoWTI$Precio_WTI
## W = 0.88542, p-value = 0.03224- Concluya sobre la validez del modelo propuesto en a)
Para concluir se puede observar que el modelo construido tiene varios inconvenientes, entre ellos que una de sus variables no se distribuye en forma normal y la otra si lo hace, generando con ello un ruido bastante alto. Por otra parte la baja respuesta del modelo a la cantidad muy reducida de 50% de los datos impide generalizar y pretender que dicho modelo explique de manera adecuada la relación entre las variables.
2. Los siguientes datos corresponden a la INFLACION y al SALARIO MINIMO LEGAL MENSUAL (SMLM) desde el año 1999 para Colombia.
La idea es establecer un modelo de regresión que ayude a determinar el comportamiento de estas dos variables tomando como variable dependiente SALARIO MINIMO LEGAL MENSUAL (SMLM) y como variable independiente INFLACION obtenga un modelo de regresión lineal simple y resuelva:
Cargar Datos
library(readr)
inflaciondf <- read_delim("D:/MAESTRIA CIENCIA DE DATOS/Segundo Semestre/Modelos Estadisticos/EntregaUnidad1/BasesdeDatos/inflacion.csv",
delim = ";", escape_double = FALSE, col_types = cols(ANO = col_number(),
INFLACION = col_number(), SMLM = col_number()),
trim_ws = TRUE)
inflaciondf## # A tibble: 17 × 3
## ANO INFLACION SMLM
## <dbl> <dbl> <dbl>
## 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.5 358000
## 7 2005 4.85 381500
## 8 2006 4.48 408000
## 9 2007 5.69 433700
## 10 2008 7.67 461500
## 11 2009 2 496900
## 12 2010 3.17 515000
## 13 2011 3.73 535600
## 14 2012 2.44 566700
## 15 2013 1.94 589500
## 16 2014 3.66 616027
## 17 2015 6.77 644350Para comenzar con el analisis realizamos un grafico de disersión junto con una li.
library(ggplot2)
ggplot(inflaciondf, aes(INFLACION,SMLM)) + geom_point()+scale_x_continuous("Inflación (%)") + scale_y_continuous("Salario mínimo legal mensual - SMLM (Pesos)") + geom_smooth()
Analisis descriptivo
attach(inflaciondf)
summary(inflaciondf)## ANO INFLACION SMLM
## Min. :1999 Min. :1.940 Min. :236460
## 1st Qu.:2003 1st Qu.:3.660 1st Qu.:332000
## Median :2007 Median :5.500 Median :433700
## Mean :2007 Mean :5.354 Mean :437079
## 3rd Qu.:2011 3rd Qu.:6.990 3rd Qu.:535600
## Max. :2015 Max. :9.230 Max. :644350hist(INFLACION,col="green")
- Escriba la ecuación del modelo de regresión lineal simple.
Una vez establecido el modelo de regresión lineal se consigue la euación:
\[SMLM=648486 - 39.49*INFLACION\]
modelo2=lm(SMLM~INFLACION,data=inflaciondf)
summary(modelo2)##
## Call:
## lm(formula = SMLM ~ INFLACION, data = inflaciondf)
##
## 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 ***
## 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.00145- Plantee y valide las hipótesis correspondientes a la linealidad general del modelo propuesto en a).
- Construya una gráfica de residuales y haga un análisis cualitativo de los supuestos del modelo propuesto en a).
Hipotesis: La inflación se distribuye normalmente. El SMLM se distribuye normalmente.
Para validar los supuestos, se realiza el analisis de residuos y se encuentra que cumple con con la media cero y es visible que no hay ninguna tendencia entre los residuales y los ajustados.
Por otra parte existe un ajuste aprobable en la recta de normalidad qqplot. Y al realizar las pruebas de Shapiro para cada variable, el p valor (0.5717) es mayor a \(\alpha=0.05\), no se rechaza que la variable INFLACION presenta un comportamiento normal. Para la variable SMLM con el el p valor 0.6218, tambien muestra un comportamiento de manera normal, confirmando la observación detallada en el analisis estadistico descriptivo de la parte inicial.
par(mfrow=c(2,2))
plot(modelo2)
shapiro.test(inflaciondf$INFLACION)##
## Shapiro-Wilk normality test
##
## data: inflaciondf$INFLACION
## W = 0.95677, p-value = 0.5717shapiro.test(inflaciondf$SMLM)##
## Shapiro-Wilk normality test
##
## data: inflaciondf$SMLM
## W = 0.95948, p-value = 0.6218- Indique e interprete el coeficiente de correlación del modelo propuesto en a).
El coeficiente de correlación que se obtiene para estas dos variables esta dado por -0.7086581, lo que ilustra una correlación negativa, mostrando que a medida que la inflación aumenta el SMLM disminuye, existiendo una corelación alta entre las dos variables y expresando que el 71% de los datos se encuentra bien caracterizado por el modelo construido.
Por otra parte tenemos un \[R^{2}=0.469\]
cor(inflaciondf$INFLACION,inflaciondf$SMLM)## [1] -0.7086581- Interprete cada uno de los coeficientes del modelo propuesto en a).
Los coeficientes encontrados en el modelo desarrollado son \(\beta_0=648486\) y \(\beta_1=-39.49\), el primero de ellos hace referencia al salario base, inicial e ideal, en el que la inflación seria 0% por otra parte el segundo obede a la tasa de reducción del salario por cada punto porcentual que suba la inflación.
- Comente sobre la conveniencia de usar el modelo propuesto en a) para predecir el SMLM para Colombia.
La utilizacion de la función \(SMLM=648486 - 39.49*INFLACION\) establece una manera de predecir cuanto sera el salario mínimo dado el porcentaje de aumento de la inflación año a año, sin embargo y aunque el modelo sea acertivo y presente buenas metricas no tiene sentido las cifras que obtenemos, es decir el SMLM no representaria en si una forma de predecir el valor de esta variable.
El módelo es mucho mas estimativo si se piensa en que lo que representa, es una perdida gradual de la perdida gradual de la capacidad de capital y consumo quepueden tener las personas que ganan un salario minimo. Existe una perdida sistematica de la capacidad adquisitiva bajo la comparación entre inflación y valor del salario minimo.
3. 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).
Cargar Datos
library(readxl)
viviendadf <- read_excel("D:/MAESTRIA CIENCIA DE DATOS/Segundo Semestre/Modelos Estadisticos/EntregaUnidad1/BasesdeDatos/Datos_Vivienda.xlsx",
col_types = c("text", "text", "numeric",
"numeric", "numeric", "numeric", "numeric",
"numeric", "text", "text", "numeric",
"numeric"))
viviendadf## # A tibble: 8,322 × 12
## Zona piso Estrato precio_…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 Zona Sur 2 6 880 237 2 5 4 Casa pance
## 2 Zona Oeste 2 4 1200 800 3 6 7 Casa miraf…
## 3 Zona Sur 3 5 250 86 NA 2 3 Apar… multi…
## 4 Zona Sur <NA> 6 1280 346 4 6 5 Apar… ciuda…
## 5 Zona Sur 2 6 1300 600 4 7 5 Casa pance
## 6 Zona Sur 3 6 513 160 2 4 4 Casa pance
## 7 Zona Sur 2 6 870 490 3 6 5 Casa ciuda…
## 8 Zona Sur 5 5 310 82.5 1 2 3 Apar… valle…
## 9 Zona Sur 9 4 240 80 1 2 3 Apar… valle…
## 10 Zona Sur 6 6 690 150 2 5 4 Apar… pance
## # … with 8,312 more rows, 2 more variables: cordenada_longitud <dbl>,
## # Cordenada_latitud <dbl>, and abbreviated variable names ¹precio_millon,
## # ²Area_contruida, ³parqueaderos, ⁴Habitaciones
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names- 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.
library(dplyr)
vivienda1df <-filter(viviendadf, Zona=="Zona Norte" & precio_millon<500 & Area_contruida<300)
head(vivienda1df,3)## # A tibble: 3 × 12
## Zona piso Estrato preci…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio corde…⁵
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl>
## 1 Zona… 2 3 135 56 1 1 3 Apar… torre… -76.5
## 2 Zona… <NA> 5 400 212 NA 2 4 Casa santa… -76.5
## 3 Zona… <NA> 3 78 54 2 1 3 Apar… chimi… -76.5
## # … with 1 more variable: Cordenada_latitud <dbl>, and abbreviated variable
## # names ¹precio_millon, ²Area_contruida, ³parqueaderos, ⁴Habitaciones,
## # ⁵cordenada_longitud
## # ℹ Use `colnames()` to see all variable namesvivienda1df## # A tibble: 1,483 × 12
## Zona piso Estrato precio_…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 Zona Norte 2 3 135 56 1 1 3 Apar… torre…
## 2 Zona Norte <NA> 5 400 212 NA 2 4 Casa santa…
## 3 Zona Norte <NA> 3 78 54 2 1 3 Apar… chimi…
## 4 Zona Norte <NA> 3 175 130 NA 3 4 Casa brisa…
## 5 Zona Norte <NA> 5 340 106 2 2 3 Apar… la fl…
## 6 Zona Norte 2 4 265 162 1 3 4 Casa zona …
## 7 Zona Norte <NA> 3 120 130 1 2 4 Casa flora…
## 8 Zona Norte <NA> 3 380 177 NA 0 0 Casa gaitan
## 9 Zona Norte 1 3 135 103 1 3 4 Apar… calim…
## 10 Zona Norte 1 3 75 54 1 2 3 Apar… calim…
## # … with 1,473 more rows, 2 more variables: cordenada_longitud <dbl>,
## # Cordenada_latitud <dbl>, and abbreviated variable names ¹precio_millon,
## # ²Area_contruida, ³parqueaderos, ⁴Habitaciones
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable namesEn este caso tenemos que hacen falta datos e información sobre los parqueaderos, en ellos es evidente que no todos las viviendas tiene la información del número exacto de parqueaderos. En este caso se encuentran 507 viviendas que no tienen información sobre parqueaderos.
sapply(vivienda1df, function(x) sum(is.na(x)))## Zona piso Estrato precio_millon
## 0 578 0 0
## Area_contruida parqueaderos Banos Habitaciones
## 0 507 0 0
## Tipo Barrio cordenada_longitud Cordenada_latitud
## 0 0 0 0Para realizar un analisis con la información completa y no poner ruido a nuestros datos eliminaremos la información de las viviendas que no entregan información sobre parqueaderos.
vivienda1df<-vivienda1df[!is.na(vivienda1df$parqueaderos),]
vivienda1df## # A tibble: 976 × 12
## Zona piso Estrato precio_…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 Zona Norte 2 3 135 56 1 1 3 Apar… torre…
## 2 Zona Norte <NA> 3 78 54 2 1 3 Apar… chimi…
## 3 Zona Norte <NA> 5 340 106 2 2 3 Apar… la fl…
## 4 Zona Norte 2 4 265 162 1 3 4 Casa zona …
## 5 Zona Norte <NA> 3 120 130 1 2 4 Casa flora…
## 6 Zona Norte 1 3 135 103 1 3 4 Apar… calim…
## 7 Zona Norte 1 3 75 54 1 2 3 Apar… calim…
## 8 Zona Norte 3 3 125 140 1 2 4 Casa calim…
## 9 Zona Norte <NA> 4 175 77 1 2 3 Apar… urban…
## 10 Zona Norte 2 3 140 130 1 2 3 Casa calim…
## # … with 966 more rows, 2 more variables: cordenada_longitud <dbl>,
## # Cordenada_latitud <dbl>, and abbreviated variable names ¹precio_millon,
## # ²Area_contruida, ³parqueaderos, ⁴Habitaciones
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names(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?).
Se observa que la distribución de los puntos en el mapa no solamente se concentra en el la zona norte de la ciudad, por el contrario existe una gran dispersión de los datos en el mapa. Toda ella debida en esencia a que aunque todos este classificados en Zona Norte las latitudes y longitudes no coinciden estrictamente con la zona norte de la ciudad de Santiago de Cali.
require(leaflet)
leaflet() %>% addCircleMarkers(lng=vivienda1df$cordenada_longitud, lat=vivienda1df$Cordenada_latitud,radius = 0.3,color = "black",label = vivienda1df$Barrio)%>% addTiles()- 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.
Analisis Exploratorio:
library(tidyverse)
library(plotly)
library(ggplot2)
grafica<- plot_ly(data=vivienda1df,x=~Area_contruida, y=~precio_millon, type="scatter", color=~Tipo)
graficagrafica1<-qplot(Estrato, precio_millon, data = vivienda1df, geom=c("boxplot"), fill=Estrato)
grafica1<-ggplotly(grafica1)
grafica1grafica2<-qplot(parqueaderos, precio_millon, data = vivienda1df, geom=c("boxplot"), fill=parqueaderos)
grafica2<-ggplotly(grafica2)
grafica2- 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 \(R^{2}\) y discuta el ajuste del modelo e implicaciones (que podrían hacer para mejorarlo).
modelovivienda <- lm(formula = precio_millon ~ Area_contruida + parqueaderos + Estrato, data = vivienda1df)
summary(modelovivienda)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + parqueaderos +
## Estrato, data = vivienda1df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -164.506 -39.027 -3.398 32.398 221.802
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -131.83094 9.66087 -13.646 <2e-16 ***
## Area_contruida 0.67929 0.03805 17.851 <2e-16 ***
## parqueaderos 28.82936 3.28318 8.781 <2e-16 ***
## Estrato 66.81995 2.25474 29.635 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 58.45 on 972 degrees of freedom
## Multiple R-squared: 0.6886, Adjusted R-squared: 0.6876
## F-statistic: 716.4 on 3 and 972 DF, p-value: < 2.2e-16Para el modelo lineal multiple establecido las tres variables son altamente significativas, mostrando que el modelo tiene una alta corelación lineal de cada una de las variables en el precio de las viviendas. El coeficiente de \(R^{2}\) y \(R^{2}\) ajustado estan muy proximos. Es muy importente observar que a medida que se aumenten variables en el modelo lineal el coeficiente \(R^{2}\) aumentará, pero se debe revisar esencialemnte el coeficiente \(R^{2}\) ajustado para saber que tanto aporta cada una de las variables, más y mas.
Para ello vamos a ir contruyendo el modelo lineal variable por variable.
modelovivienda1 <- lm(formula = precio_millon ~ Area_contruida, data = vivienda1df)
summary(modelovivienda1)##
## Call:
## lm(formula = precio_millon ~ Area_contruida, data = vivienda1df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -223.878 -72.368 -4.719 59.722 220.912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 152.38674 6.11672 24.91 <2e-16 ***
## Area_contruida 1.09067 0.04906 22.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 85.22 on 974 degrees of freedom
## Multiple R-squared: 0.3366, Adjusted R-squared: 0.336
## F-statistic: 494.3 on 1 and 974 DF, p-value: < 2.2e-16modelovivienda2 <- lm(formula = precio_millon ~ Area_contruida + Estrato, data = vivienda1df)
summary(modelovivienda2)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato, data = vivienda1df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -199.260 -40.554 -4.884 31.684 220.043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -124.43295 9.99336 -12.45 <2e-16 ***
## Area_contruida 0.81592 0.03606 22.63 <2e-16 ***
## Estrato 70.67870 2.29636 30.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.69 on 973 degrees of freedom
## Multiple R-squared: 0.6639, Adjusted R-squared: 0.6632
## F-statistic: 960.9 on 2 and 973 DF, p-value: < 2.2e-16Al realizarlo variable a variable es evidente que cada variable va aumentando el \(R^{2}\) ajustado y se evidencia que iempre hay valores grandes de significancia entre las variables seleccionadas con el precio de la vivienda.
- 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).
Se realiza una validaciión de supuestos sobre el caso de la regresión lineal completa con respecto a las tres variables.
library(ggplot2)
library(gridExtra)
plot1 <- ggplot(data = vivienda1df, aes(Area_contruida, modelovivienda$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = vivienda1df, aes(parqueaderos, modelovivienda$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = vivienda1df, aes(Estrato, modelovivienda$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3)
qqnorm(modelovivienda$residuals)
qqline(modelovivienda$residuals)
shapiro.test(modelovivienda$residuals)##
## Shapiro-Wilk normality test
##
## data: modelovivienda$residuals
## W = 0.99025, p-value = 4.54e-06Tanto el análisis gráfico como es test de hipótesis confirman la normalidad. Y se tiene una variabilidad constante de los residuos lo que se denomina homocedasticidad. Al representar los residuos frente a los valores ajustados por el modelo, los primeros se tienen que distribuir de forma aleatoria en torno a cero, manteniendo aproximadamente la misma variabilidad a lo largo del eje X. S No se observan patrones especificos.
Sin embargo todo el analisis propone variables númericas. Pero las variables parqueaderos y estrato son realmente categoricas.
Se podria mejorar el modelo tratando de construir un modelo lineal que contemple variables categoricas.
- 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?
preciovivienda=-131.83094+ 0.67929 *(100)+ 28.82936*0 + 66.81995*4
preciovivienda## [1] 203.3779preciovivienda1=-131.83094+ 0.67929 *(100)+ 28.82936*1 + 66.81995*4
preciovivienda1## [1] 232.2072preciovivienda2=-131.83094+ 0.67929 *(100)+ 28.82936*2 + 66.81995*4
preciovivienda2## [1] 261.0366preciovivienda3=-131.83094+ 0.67929 *(100)+ 28.82936*3 + 66.81995*4
preciovivienda3## [1] 289.8659preciovivienda4=-131.83094+ 0.67929 *(100)+ 28.82936*4 + 66.81995*4
preciovivienda4## [1] 318.6953preciovivienda5=-131.83094+ 0.67929 *(100)+ 28.82936*5 + 66.81995*4
preciovivienda5## [1] 347.5247preciovivienda6=-131.83094+ 0.67929 *(100)+ 28.82936*6 + 66.81995*4
preciovivienda6## [1] 376.354En nuestro modelo la cantidad del parquedero afecta considerablemente el costo de la vivienda es por ello que para el ejemplo propuesto si la vivienda tiene bastantes parqueaderos puede ser una vivienda en una buena oferta, sin embargo para el modelo construido no es totalmente buena la oferta.
Si utilizamos el modelo que implica solo la relación entre precio de vivienda en función de area construida y estrato tambien observamos que la vivienda ofrecida es bastante costosa. No es una buena oferta.
preciovivienda= -124.43295 + 0.81592 * 100 + 70.67870*4
preciovivienda## [1] 239.8739- 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.
Para realizar este analisis podemos realizar un filtro de viviendas con estas carateristicas que no superen los 150 millones.
library(dplyr)
vivienda2df <-filter(vivienda1df, Zona=="Zona Norte" & Estrato==4 & Area_contruida>100 & precio_millon <400)
head(vivienda2df,5)## # A tibble: 5 × 12
## Zona piso Estrato preci…¹ Area_…² parqu…³ Banos Habit…⁴ Tipo Barrio corde…⁵
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl>
## 1 Zona… 2 4 265 162 1 3 4 Casa zona … -76.5
## 2 Zona… <NA> 4 250 115 1 3 3 Casa villa… -76.5
## 3 Zona… 2 4 330 165 1 4 4 Casa el bo… -76.5
## 4 Zona… 2 4 370 240 2 5 5 Casa el bo… -76.5
## 5 Zona… 2 4 350 280 2 3 4 Casa la me… -76.5
## # … with 1 more variable: Cordenada_latitud <dbl>, and abbreviated variable
## # names ¹precio_millon, ²Area_contruida, ³parqueaderos, ⁴Habitaciones,
## # ⁵cordenada_longitud
## # ℹ Use `colnames()` to see all variable namesrequire(leaflet)
leaflet() %>% addCircleMarkers(lng=vivienda2df$cordenada_longitud, lat=vivienda2df$Cordenada_latitud,radius = 0.3,color = "black",label = vivienda1df$Barrio)%>% addTiles()- 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.
Cargar Datos
library(readxl)
arbolesdf <- read_excel("D:/MAESTRIA CIENCIA DE DATOS/Segundo Semestre/Modelos Estadisticos/EntregaUnidad1/BasesdeDatos/data_arboles.xlsx",
col_types = c("text", "text", "numeric",
"numeric", "numeric"))
arbolesdf## # A tibble: 90 × 5
## finca mg peso diametro altura
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 FINCA_1 GENOTIPO_1 13.7 4.7 5
## 2 FINCA_1 GENOTIPO_1 14.6 5.3 5.6
## 3 FINCA_1 GENOTIPO_1 15.9 4.8 5.8
## 4 FINCA_1 GENOTIPO_1 8.99 3.2 4.3
## 5 FINCA_1 GENOTIPO_1 6.99 2.2 3.3
## 6 FINCA_1 GENOTIPO_2 19.3 6.3 7.9
## 7 FINCA_1 GENOTIPO_2 21.4 6.6 8.3
## 8 FINCA_1 GENOTIPO_2 13.8 5.3 7.3
## 9 FINCA_1 GENOTIPO_2 11.9 4.9 6.7
## 10 FINCA_1 GENOTIPO_2 16.6 5.9 7.1
## # … with 80 more rows
## # ℹ Use `print(n = ...)` to see more rowsSe puede buscar esencialmente corelaciones graficas entre las diferentes variables, y en ellas se observa que el peso puede estar correlacionado con el diametro y la altura en algun grado leve o mederado.
#arbolesdf$peso <- as.factor(arbolesdf$peso)
#pairs(x = arbolesdf)Si hacemos el caculo de correlacion entre las variables conseguimos una correlación alta entre ellas.
cor(arbolesdf$peso,arbolesdf$diametro)## [1] 0.908123cor(arbolesdf$peso,arbolesdf$altura)## [1] 0.8582009Adicionalmente podemos evaluar la significancia de la corelación.
cor.test(arbolesdf$peso,arbolesdf$diametro)##
## Pearson's product-moment correlation
##
## data: arbolesdf$peso and arbolesdf$diametro
## t = 20.346, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8634081 0.9386817
## sample estimates:
## cor
## 0.908123cor.test(arbolesdf$peso,arbolesdf$altura)##
## Pearson's product-moment correlation
##
## data: arbolesdf$peso and arbolesdf$altura
## t = 15.684, df = 88, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7918402 0.9045332
## sample estimates:
## cor
## 0.8582009En ambos casos se observa que hay un valor de significancia bastante relevante y que en el intervalo de confianza del 95% ambas variables explican desde un 80% a un 94% de los datos.
Ajustando el modelo de regresión multiple se consigue:
modelo5 <- lm(peso ~ diametro + altura, data = arbolesdf)
summary(modelo5)##
## Call:
## lm(formula = peso ~ diametro + altura, data = arbolesdf)
##
## 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-16Se consigue dentro del modelo un valor de significancia bastante alto para el diametro, sin embargo no lo es tan bueno para la corelación con la altura. Adicionalmente conseguimos un buen \(R^{2}\), mostrando una alta efectividad en el modelo.
Modelo con los datos completos.
\[peso=-9.1205 + 4.7395*diametro + 0.3132 *altura\]
Validación Simple
Se selecciona un subconjunto de dos datos para hacer una validación simple.
library(ISLR)
data(arbolesdf)
dim(arbolesdf)## [1] 90 5set.seed(1)
train <- sample(x = 1:90, 45)
modelo6 <- lm(peso ~ diametro + altura, data = arbolesdf, subset=train)
summary(modelo6)##
## Call:
## lm(formula = peso ~ diametro + altura, data = arbolesdf, subset = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.4847 -2.0586 0.5223 1.6020 10.3523
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.3593 1.9115 -4.896 1.49e-05 ***
## diametro 3.4659 0.9974 3.475 0.0012 **
## altura 1.4493 0.7642 1.896 0.0648 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.258 on 42 degrees of freedom
## Multiple R-squared: 0.8631, Adjusted R-squared: 0.8566
## F-statistic: 132.4 on 2 and 42 DF, p-value: < 2.2e-16predicciones <-predict(object = modelo6, newdata = arbolesdf[-train, ])
predicciones## 1 2 3 4 5 6 7 8
## 15.682948 7.963543 3.048336 19.589823 17.333886 21.379531 4.667825 15.191427
## 9 10 11 12 13 14 15 16
## 3.338190 10.131293 24.183728 24.416843 24.215176 17.768667 9.784699 17.245698
## 17 18 19 20 21 22 23 24
## 14.674615 20.427937 12.185563 22.482208 9.841439 25.979593 20.950906 20.131926
## 25 26 27 28 29 30 31 32
## 17.069323 23.314165 23.547280 21.001489 31.966029 34.707330 15.248167 15.336354
## 33 34 35 36 37 38 39 40
## 24.788728 22.249093 22.015978 27.655821 25.891405 21.959239 15.046500 17.705771
## 41 42 43 44 45
## 12.009188 17.415916 11.517667 15.506573 13.628678error<-mean((arbolesdf$peso[-train] - predicciones)^2)
error## [1] 14.39137Obteniendo un error MSE de 14.39137.
Por otra parte se realiza un proceso de validación cruzada con 10 grupos.
Validación Cruzada
library(boot)
# Se genera el modelo lineal
modelo7 <- glm(peso ~ diametro + altura, data = arbolesdf)
# Se emplea la función cv.glm() para la validación, empleando en este caso k=10
set.seed(1)
cv_error <- cv.glm(data = arbolesdf, glmfit = modelo7, K = 10)
cv_error$delta## [1] 12.75600 12.68806predicciones1 <-predict(object = modelo7, newdata = arbolesdf[-train, ])
predicciones1## 1 2 3 4 5 6 7 8
## 15.445213 7.392339 2.339721 18.284683 16.201003 21.065725 3.886823 14.939951
## 9 10 11 12 13 14 15 16
## 2.402353 9.856016 24.065980 23.717299 23.274669 16.294952 9.382071 16.581001
## 17 18 19 20 21 22 23 24
## 13.232065 19.263890 11.497066 20.904940 9.793384 24.853087 18.977841 21.195194
## 25 26 27 28 29 30 31 32
## 17.340996 23.878083 23.529402 21.383090 30.536231 35.119108 15.351264 14.971267
## 33 34 35 36 37 38 39 40
## 25.393869 21.253622 21.602303 26.811502 25.233084 21.190989 14.908635 17.877573
## 41 42 43 44 45
## 12.257061 17.814941 11.751799 16.205207 13.804163Obteniendo un error RMSE de:
Metrica RMSE
Para la metrica RMSE que se define como:
\[RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i-<y_i>)^{2}}\]
RMSE<-sqrt(mean((arbolesdf$peso[-train] - predicciones1)^2))
RMSE## [1] 3.52681y un error MAE de:
Metrica MAE
Para la metrica MAE que se define como:
\[MAE=\frac{1}{N}\sum_{i=1}^{N}|y_i-<y_i>|\]
MAE<-mean(abs(arbolesdf$peso[-train] - predicciones1))
MAE## [1] 2.92943