Informe de aplicación de regresión lineal

https://rpubs.com/sgcifuentes/935217

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

Proponga un modelo de regresión lineal simple que permita predecir el valor de las Acciones de Ecopetrol con base en el Precio del barril de petróleo en Colombia. Indique la ecuación de regresión y el valor del R2.

require(ggplot2)
require(plotly)

prec_ecop=c(1090, 1170, 1160, 1230, 1155, 1165, 1205, 1170, 1150, 1130, 1110, 1105, 
1085, 1060, 1035, 1015, 955, 961) 

prec_petr=c(35.62,36.31,37.35,34.95,34.53,35.81,36.14,37.50,37.80,36.81,37.87,37.04,
36.76,35.97,33.97,33.27,31.41,30.44)

bd=data.frame(prec_ecop,prec_petr)

#Realizamos analisis exploratorio bivariado -(Relación entre precio acciones ecopetrol contra precio del petroleo WTI)

Graf_tendencia_accion_Ecopetrol = ggplot(data = bd, mapping = aes(x=prec_petr, y=prec_ecop)) + geom_point()+theme_bw()+geom_smooth()
ggplotly(Graf_tendencia_accion_Ecopetrol)

cor(bd$prec_petr,bd$prec_ecop)

## [1] 0.7074373

Se observa que la relación entre el precio del petroleo WTI y el precio de las acciones de Ecopetrol es fuerte; Se evidencia en el coeficiente de correlación de 0.7, es decir, que el incremento del WTI tiene influencia directa en el precio de la acción de ecopetrol.

ESTIMACIÓN DE MODELO DE REGRESIÓN LINEAL SIMPLE

modelo_ecopetrol=lm(prec_ecop ~ prec_petr,data=bd)
modelo_ecopetrol

## 
## Call:
## lm(formula = prec_ecop ~ prec_petr, data = bd)
## 
## Coefficients:
## (Intercept)    prec_petr  
##      177.77        26.19

El modelo estimado es: Precio_Acción_Ecopetrol = 177.77 + (26.19 * precio_petroleo_WTI), donde Bo=177.77 y B1=26.19. De esto se puede interpretar que si el precio de la acción de ecopetrol no dependiera un 70% del precio del petroleo WTI, tendría un costo de 177.77 pesos, sin embargo, este precio de la acción varía 26.19 veces lo que cueste el precio del petroleo WTI.

# Realizamos el Summary del modelo

summary(modelo_ecopetrol)

## 
## Call:
## lm(formula = prec_ecop ~ prec_petr, data = bd)
## 
## 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   
## prec_petr     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.001024

Se observa en la tabla Summary del modelo que el precio del petroleo WTI es significativo en el modelo porque el valor p lo indica con un nivel de confianza del 99.8% y además el modelo logra explicar el 50% de la variablidad del precio de la acción de Ecopetrol con el R2 de 0.5005

b.Pruebe la significancia del modelo propuesto en “a)” plantee las hipótesis respectivas y use el concepto de Valor _p para tomar la decisión sobre las hipótesis. Use α = 0.05

Las Hipotesis se expresarian de la siguiente manera: Hipotesis nula (Ho): B1 = 0 (no hay relación) Hipotesis alterna (H1): B1 ≠ 0 (Si hay relación)

Dado que se tiene un p-value: 0.001024 en el modelo y este es inferior a 0.05, se rechaza la hipotesis nula y se concluye lo siguiente: Si existe una relación significativa entre la variable dependiente que es el precio de la acción del petroleo y la variable independiente que es el precio del petroleo WTI.

El modelo estimado es: Precio_Acción_Ecopetrol = 177.77 + (26.19 * precio_petroleo_WTI), donde Bo=177.77 y B1=26.19. De esto se puede interpretar que si el precio de la acción de ecopetrol no dependiera un 70% del precio del petroleo WTI, tendría un costo de 177.77 pesos, sin embargo, este precio de la acción varía 26.19 veces lo que cueste el precio del petroleo WTI.

Análsis de residuos

Se puede observar en la gráfica de normalidad de los residuos que estos están bien aproximados a la recta, es decir, hay una mínima cantidad de error en el modelo y los residuos se consideran que tienen normalidad.

residuos = residuals(modelo_ecopetrol)
plot(prec_petr,residuos);abline(h=0)

qqnorm(residuos);qqline(residuos)

res=modelo_ecopetrol$residuals

mean(res)

## [1] -5.527407e-15

t.test(res)

## 
##  One Sample t-test
## 
## data:  res
## t = -4.2309e-16, df = 17, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -27.56364  27.56364
## sample estimates:
##     mean of x 
## -5.527407e-15

shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.89259, p-value = 0.04276

lmtest::gqtest(modelo_ecopetrol)

## 
##  Goldfeld-Quandt test
## 
## data:  modelo_ecopetrol
## GQ = 0.17924, df1 = 7, df2 = 7, p-value = 0.9813
## alternative hypothesis: variance increases from segment 1 to 2

lmtest::dwtest(modelo_ecopetrol)

## 
##  Durbin-Watson test
## 
## data:  modelo_ecopetrol
## DW = 0.74504, p-value = 0.0004666
## alternative hypothesis: true autocorrelation is greater than 0

En el test Simple se valida que el valor p es 1, es decir, Ok En el test de Shapiro el valor p da 0.04, es decir, correcto se rechaza la hipotesis nula. En el test Goldfeld el valor de p es de 0.98 es decir que las varianzas se asumen iguales. En el test de Durbin Watson el valor de p es 0.00046, es decir que no hay independencia en los errores.

Se validan todos los supuestos y se concluye que todos estos se cumplen y de esta forma el modelo está listo para cualquier uso de inferencias estadísticas.

El modelo estimado es valido y cumple con todas las condiciones para ser usado de acuerdo a lo siguiente: Precio_Acción_Ecopetrol = 177.77 + (26.19 * precio_petroleo_WTI), donde Bo=177.77 y B1=26.19. De esto se puede interpretar que si el precio de la acción de ecopetrol no dependiera un 70% del precio del petroleo WTI, tendría un costo de 177.77 pesos, sin embargo, este precio de la acción varía 26.19 veces lo que cueste el precio del petroleo WTI.

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:

Escriba la ecuación del modelo de regresión lineal simple

infl=c(9.23, 8.75, 7.65, 6.99, 6.49, 5.50, 4.85, 4.48, 5.69, 7.67, 2.00, 3.17, 
3.73, 2.44,  1.94, 3.66, 6.77)  
  
smlm=c(236460, 260100, 286000, 309000, 332000, 358000, 381500, 408000, 433700, 
461500, 496900, 515000, 535600, 566700, 589500, 616027, 644350)

bdSML=data.frame(infl,smlm)
head(bdSML)

##   infl   smlm
## 1 9.23 236460
## 2 8.75 260100
## 3 7.65 286000
## 4 6.99 309000
## 5 6.49 332000
## 6 5.50 358000

cor(bdSML$infl,bdSML$smlm)

## [1] -0.7086581

modelo_salario_minimo=lm(smlm ~ infl,data=bdSML)
modelo_salario_minimo

## 
## Call:
## lm(formula = smlm ~ infl, data = bdSML)
## 
## Coefficients:
## (Intercept)         infl  
##      648486       -39489

summary(modelo_salario_minimo)

## 
## Call:
## lm(formula = smlm ~ infl, data = bdSML)
## 
## 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 ***
## infl          -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

Salario Minimo = 648.486 + (-39489 * inflacion), donde Bo=648.486 y B1= -39489

Esto indica que si no hubiera inflación el salario minimo sería 648.486 y que por cada punto de aumento en inflación el salario minimo disminuye 39.489

Las Hipotesis se expresarian de la siguiente manera: Hipotesis nula (Ho): B1 = 0 (no hay relación) Hipotesis alterna (H1): B1 ≠ 0 (Si hay relación)
Se observa que la relación entre la inflación y el salario minimo es fuerte pero a la inversa; El coeficiente de correlación es -0.7, esto indica hay una relación inversa directa entre ambas variables y que por cada punto de aumento en inflación el salario minimo disminuye 39.489
El modelo estimado es Salario Minimo = 648.486 + (-39489 * inflacion), donde Bo=648.486 y B1= -39489. De esto se puede interpretar que si ho existiera una relación del 50% entre ambas variables, el salario minimo no se vería afectado por la inflación, pero tenemos un valor que nos confirma dicha hipotesis que es el valor p.

Se observa en la tabla Summary del modelo que la inflación es significativo en el modelo porque el valor p lo indica con un nivel de confianza del 99.8% y además el modelo logra explicar el 50% de la variablidad del salario minimo con el R2 de 0.5022.

residuos_SML = residuals(modelo_salario_minimo)
plot(infl,residuos_SML);abline(h=0)

qqnorm(residuos_SML);qqline(residuos_SML)

res_SML=modelo_salario_minimo$residuals

mean(res_SML)

## [1] -1.491304e-12

t.test(res_SML)

## 
##  One Sample t-test
## 
## data:  res_SML
## t = -6.7462e-17, df = 16, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -46862.45  46862.45
## sample estimates:
##     mean of x 
## -1.491304e-12

shapiro.test(res_SML)

## 
##  Shapiro-Wilk normality test
## 
## data:  res_SML
## W = 0.78826, p-value = 0.001407

lmtest::gqtest(modelo_salario_minimo)

## 
##  Goldfeld-Quandt test
## 
## data:  modelo_salario_minimo
## GQ = 140.68, df1 = 7, df2 = 6, p-value = 3.171e-06
## alternative hypothesis: variance increases from segment 1 to 2

lmtest::dwtest(modelo_salario_minimo)

## 
##  Durbin-Watson test
## 
## data:  modelo_salario_minimo
## DW = 0.68432, p-value = 0.0002714
## alternative hypothesis: true autocorrelation is greater than 0

En el test Simple se valida que el valor p es 1, es decir, Ok En el test de Shapiro el valor p da 0.001, es decir, correcto se rechaza la hipotesis nula. En el test Goldfeld el valor de p es de 3.171e-06 es decir que las varianzas no son iguales, sin embargo no afecta el modelo. En el test de Durbin Watson el valor de p es 0.00027, es decir que no hay independencia en los errores.

Se validan todos los supuestos y se concluye que todos estos se cumplen y de esta forma el modelo está listo para cualquier uso de inferencias estadísticas.

predict (modelo_salario_minimo,list(infl=8.75),interval = "confidence",level =0.95)

##        fit    lwr      upr
## 1 302954.3 214813 391095.6

El modelo es conveniente y validando la data aleatoriamente se observa que los datos están dentro del rango predecido.

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 .

Realice un filtro a la base de datos e incluya solo las ofertas de apartamentos, de la zona norte de la ciudad con precios inferiores a los 500 millones de pesos y áreas menores a 300 mt2. Presente los primeros 3 registros de la base y algunas tablas que comprueben la consulta. (Adicional un mapa con los puntos de la base, discutir si todos los puntos se ubican en la zona norte o se presentan valores en otras zonas, por que?).

A continuación podemos observar la base de datos filtrada por apartamentos de la zona norte de la ciudad con precios inferiores a los 500 millones de pesos y áreas menores a 300 mt2 y al graficar estos lugares nos encontramos que no todos los puntos se concentran en la zona norte y esto se debe a que la longitud y latitud está errada, es decir, el vendedor puso mal la dirección en la plataforma.

library(readxl)
datos = read_excel("C:/Users/User/Downloads/Datos_Vivienda (1).xlsx")

#View(datos)
ID=1:dim(datos)[1]

datos_Id=data.frame(ID,datos)

#View(datos_Id)

require(dplyr)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

db_filtrada2 = datos_Id %>% filter(Tipo == "Apartamento" & Zona =="Zona Norte" & precio_millon < "500" & Area_contruida < "300")

#View(db_filtrada)

db_filtrada2[1:3,]

##    ID       Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1  89 Zona Norte   NA       5           340            106            2     2
## 2 104 Zona Norte    1       3           135            103            1     3
## 3 183 Zona Norte   NA       5           390            102            2     3
##   Habitaciones        Tipo                 Barrio cordenada_longitud
## 1            3 Apartamento               la flora          -76.48200
## 2            4 Apartamento          calimio norte          -76.48347
## 3            3 Apartamento urbanizaciv=n la flora          -76.48800
##   Cordenada_latitud
## 1           3.43500
## 2           3.48626
## 3           3.46400

require(leaflet)

## Loading required package: leaflet

## Warning: package 'leaflet' was built under R version 4.1.3

leaflet() %>% addCircleMarkers(lng = db_filtrada2$cordenada_longitud,lat = db_filtrada2$Cordenada_latitud,radius = 0.3,color = "black") %>% 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.

precio_promedio=mean(db_filtrada2$precio_millon,na.rm = TRUE)
mediana_precio=median(db_filtrada2$precio_millon,na.rm = TRUE)
promedio_area=mean(db_filtrada2$Area_contruida,na.rm = TRUE)
cantidad_ofertas=length(db_filtrada2$Area_contruida)

Resultado = data.frame(precio_promedio, mediana_precio, promedio_area, cantidad_ofertas)
Resultado

##   precio_promedio mediana_precio promedio_area cantidad_ofertas
## 1        370.5778            350      131.8952              315

require(ggplot2)
require(plotly)
attach(db_filtrada2)

## The following object is masked _by_ .GlobalEnv:
## 
##     ID

Graf_area_construida = ggplot(data = db_filtrada2, mapping = aes(x=Area_contruida, y=precio_millon)) + geom_point()+theme_bw()+geom_smooth()

ggplotly(Graf_area_construida)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

cor(db_filtrada2$Area_contruida,db_filtrada2$precio_millon)

## [1] 0.4764374

Se observa que la relación entre el área construida y el precio de los apartamentos es media; Se evidencia en el coeficiente de correlación de 0.47, es decir, el área construida si tiene influencia directa en el precio de la vivienda.

Graf_estrato = ggplot(data = db_filtrada2, mapping = aes(x=Estrato, y=precio_millon)) + geom_point()+theme_bw()+geom_smooth()

ggplotly(Graf_estrato)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 2.985

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.015

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 4.9425e-016

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1.0302

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 2.985

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.015

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 4.9425e-016

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1.0302

cor(db_filtrada2$Estrato,db_filtrada2$precio_millon)

## [1] 0.4657815

db_filtrada2$parqueaderos = as.numeric(db_filtrada2$parqueaderos)

## Warning: NAs introducidos por coerción

db_filtrada2$parqueaderos[is.na(db_filtrada2$parqueaderos)] = 0

View(db_filtrada2)

Graf_parqueadero = ggplot(data = db_filtrada2, mapping = aes(x=parqueaderos, y=precio_millon)) + geom_point()+theme_bw()+geom_smooth()

ggplotly(Graf_parqueadero)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.02

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.02

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.1342e-015

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## -0.02

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2.02

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 1.1342e-015

## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1

cor(db_filtrada2$parqueaderos,db_filtrada2$precio_millon)

## [1] 0.1887774

Se observa que la relación entre el estrato y el precio de los apartamentos es media; Se evidencia en el coeficiente de correlación de 0.46, es decir, el estrato si tiene influencia directa en el precio de la vivienda, sin embargo, los parqueaderos que tenga la vivienda no influencian en el precio de la vivienda, esto se refleja en el coeficiente de correlación de 0.18.

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).

cor(db_filtrada2[4:7])

##                  Estrato precio_millon Area_contruida parqueaderos
## Estrato        1.0000000     0.4657815     0.14831835   0.18149309
## precio_millon  0.4657815     1.0000000     0.47643741   0.18877736
## Area_contruida 0.1483184     0.4764374     1.00000000  -0.01609626
## parqueaderos   0.1814931     0.1887774    -0.01609626   1.00000000

plot(x = db_filtrada2[4:7])

modelo_vivienda_1 = lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, data = db_filtrada2)
modelo_vivienda_1

## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = db_filtrada2)
## 
## Coefficients:
##    (Intercept)  Area_contruida         Estrato    parqueaderos  
##       -427.410           1.794         105.047          26.427

step(modelo_vivienda_1, direction = "both",trace = 0)

## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = db_filtrada2)
## 
## Coefficients:
##    (Intercept)  Area_contruida         Estrato    parqueaderos  
##       -427.410           1.794         105.047          26.427

summary(modelo_vivienda_1)

## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = db_filtrada2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -358.74  -75.77   -3.69   57.64  686.21 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -427.4097    62.9297  -6.792 5.63e-11 ***
## Area_contruida    1.7944     0.1887   9.510  < 2e-16 ***
## Estrato         105.0467    12.4673   8.426 1.35e-15 ***
## parqueaderos     26.4269     9.3189   2.836  0.00487 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 124.1 on 311 degrees of freedom
## Multiple R-squared:  0.4021, Adjusted R-squared:  0.3963 
## F-statistic: 69.71 on 3 and 311 DF,  p-value: < 2.2e-16

plot(modelo_vivienda_1)

data=cbind(db_filtrada2[5], db_filtrada2[,4:7])
Mcor=cor(data)
Mcor

##                precio_millon   Estrato precio_millon Area_contruida
## precio_millon      1.0000000 0.4657815     1.0000000     0.47643741
## Estrato            0.4657815 1.0000000     0.4657815     0.14831835
## precio_millon      1.0000000 0.4657815     1.0000000     0.47643741
## Area_contruida     0.4764374 0.1483184     0.4764374     1.00000000
## parqueaderos       0.1887774 0.1814931     0.1887774    -0.01609626
##                parqueaderos
## precio_millon    0.18877736
## Estrato          0.18149309
## precio_millon    0.18877736
## Area_contruida  -0.01609626
## parqueaderos     1.00000000

corrplot::corrplot(Mcor, method = "number")

#library(corrplot)

#corrplot.mixed(corr = cor(db_filtrada2[,c("Area_construida","Estrato","parqueaderos")] method = "pearson"))

Se observa que las variables que más impacto tienen en el precio del apartamento son Area Construida y Estrato con una correlación de 0.47 y no existe colinealidad entre las demás variables.

Una Estimación adicional

db_filtrada2$lnprecio_millon = log(db_filtrada2$precio_millon)
db_filtrada2$lnArea_contruida = log(db_filtrada2$Area_contruida)

ggplot(data = db_filtrada2, aes(x=lnArea_contruida, y=lnprecio_millon) )+ 
           geom_point() + 
           theme_bw() +
           geom_smooth(method='lm', formula=y~x, se=FALSE, col='dodgerblue1')

modelo_2 = lm(log(db_filtrada2$precio_millon) ~ log(db_filtrada2$Area_contruida), data = db_filtrada2)
summary(modelo_2)

## 
## Call:
## lm(formula = log(db_filtrada2$precio_millon) ~ log(db_filtrada2$Area_contruida), 
##     data = db_filtrada2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06355 -0.19745 -0.00487  0.21294  1.01221 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       3.32582    0.34696   9.586  < 2e-16 ***
## log(db_filtrada2$Area_contruida)  0.52121    0.07146   7.294 2.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3075 on 313 degrees of freedom
## Multiple R-squared:  0.1453, Adjusted R-squared:  0.1425 
## F-statistic:  53.2 on 1 and 313 DF,  p-value: 2.485e-12

Comparación entre modelo 1 y 2

library(memisc)

## Warning: package 'memisc' was built under R version 4.1.3

## Loading required package: lattice

## Loading required package: MASS

## Warning: package 'MASS' was built under R version 4.1.3

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:plotly':
## 
##     select

## 
## Attaching package: 'memisc'

## The following objects are masked from 'package:dplyr':
## 
##     collect, recode, rename, syms

## The following objects are masked from 'package:plotly':
## 
##     rename, style

## The following object is masked from 'package:ggplot2':
## 
##     syms

## The following objects are masked from 'package:stats':
## 
##     contr.sum, contr.treatment, contrasts

## The following object is masked from 'package:base':
## 
##     as.array

memisc::mtable(modelo_vivienda_1, modelo_2)

## 
## Calls:
## modelo_vivienda_1: lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = db_filtrada2)
## modelo_2: lm(formula = log(db_filtrada2$precio_millon) ~ log(db_filtrada2$Area_contruida), 
##     data = db_filtrada2)
## 
## ========================================================================================
##                                     modelo_vivienda_1             modelo_2              
##                                       -------------    -------------------------------  
##                                       precio_millon    log(db_filtrada2$precio_millon)  
## ----------------------------------------------------------------------------------------
##   (Intercept)                          -427.410***                 3.326***             
##                                         (62.930)                  (0.347)               
##   Area_contruida                          1.794***                                      
##                                          (0.189)                                        
##   Estrato                               105.047***                                      
##                                         (12.467)                                        
##   parqueaderos                           26.427**                                       
##                                          (9.319)                                        
##   log(db_filtrada2$Area_contruida)                                 0.521***             
##                                                                   (0.071)               
## ----------------------------------------------------------------------------------------
##   R-squared                               0.402                    0.145                
##   N                                     315                      315                    
## ========================================================================================
##   Significance: *** = p < 0.001; ** = p < 0.01; * = p < 0.05

modelo_vivienda_3 = lm(formula = precio_millon ~ Area_contruida + Estrato, data = db_filtrada2)
modelo_vivienda_3

## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato, data = db_filtrada2)
## 
## Coefficients:
##    (Intercept)  Area_contruida         Estrato  
##       -420.018           1.771         111.621

summary(modelo_vivienda_3)

## 
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato, data = db_filtrada2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -366.19  -73.44   -9.33   64.06  644.79 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -420.0176    63.5813  -6.606  1.7e-10 ***
## Area_contruida    1.7707     0.1906   9.290  < 2e-16 ***
## Estrato         111.6214    12.3873   9.011  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.5 on 312 degrees of freedom
## Multiple R-squared:  0.3866, Adjusted R-squared:  0.3827 
## F-statistic: 98.33 on 2 and 312 DF,  p-value: < 2.2e-16

memisc::mtable(modelo_vivienda_1, modelo_2, modelo_vivienda_3)

## 
## Calls:
## modelo_vivienda_1: lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos, 
##     data = db_filtrada2)
## modelo_2: lm(formula = log(db_filtrada2$precio_millon) ~ log(db_filtrada2$Area_contruida), 
##     data = db_filtrada2)
## modelo_vivienda_3: lm(formula = precio_millon ~ Area_contruida + Estrato, data = db_filtrada2)
## 
## ===========================================================================================================
##                                     modelo_vivienda_1             modelo_2              modelo_vivienda_3  
##                                       -------------    -------------------------------    -------------    
##                                       precio_millon    log(db_filtrada2$precio_millon)    precio_millon    
## -----------------------------------------------------------------------------------------------------------
##   (Intercept)                          -427.410***                 3.326***                -420.018***     
##                                         (62.930)                  (0.347)                   (63.581)       
##   Area_contruida                          1.794***                                            1.771***     
##                                          (0.189)                                             (0.191)       
##   Estrato                               105.047***                                          111.621***     
##                                         (12.467)                                            (12.387)       
##   parqueaderos                           26.427**                                                          
##                                          (9.319)                                                           
##   log(db_filtrada2$Area_contruida)                                 0.521***                                
##                                                                   (0.071)                                  
## -----------------------------------------------------------------------------------------------------------
##   R-squared                               0.402                    0.145                      0.387        
##   N                                     315                      315                        315            
## ===========================================================================================================
##   Significance: *** = p < 0.001; ** = p < 0.01; * = p < 0.05

De acuerdo a los resultados de los 3 modelos, se concluye que el modelo 1 representa en un 40% la razón por la cual el precio de un apartamento tiene mayor precio de acuerdo a su Area construida, estrato y parqueadero que tenga. La variable parqueader a pesar de ser la menos relevante, genera una variación entre el modelo y el modelo 3 de 2%.

t.test(modelo_vivienda_1$residuals)

## 
##  One Sample t-test
## 
## data:  modelo_vivienda_1$residuals
## t = -2.9355e-16, df = 314, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -13.69376  13.69376
## sample estimates:
##     mean of x 
## -2.043053e-15

shapiro.test(modelo_vivienda_1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_vivienda_1$residuals
## W = 0.87355, p-value = 2.06e-15

lmtest::bptest(modelo_vivienda_1)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_vivienda_1
## BP = 110.33, df = 3, p-value < 2.2e-16

lmtest::dwtest(modelo_vivienda_1)

## 
##  Durbin-Watson test
## 
## data:  modelo_vivienda_1
## DW = 1.7168, p-value = 0.005548
## alternative hypothesis: true autocorrelation is greater than 0

En el test Simple se valida que el valor p es 1, es decir, Ok

En el test de Shapiro el valor p da 0.001, es decir, correcto se rechaza la hipotesis nula.

En el test de Durbin Watson el valor de p es 0.00027, es decir que no hay independencia en los errores.

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?

View(db_filtrada2)

predict(modelo_vivienda_1, list(Area_contruida=100,Estrato=4,parqueaderos=1),interval = "confidence", level = 0.95)

##        fit      lwr      upr
## 1 198.6441 169.8866 227.4016

De acuerdo al análisis del modelo el precio de este apartamento tiene un rango de precio entre 169 y 227 millones, por lo que la oferta de 450 millones es una oferta muy costosa, le haría una oferta de 160 millones al vendedor.

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.

require(dplyr)

db_filtrada3 = datos_Id %>% filter(Tipo == "Apartamento" & Zona =="Zona Norte" & precio_millon < "400" & Area_contruida == "100" & parqueaderos == "1")

View(db_filtrada3)

db_filtrada3[1:5,]

##     ID       Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1  616 Zona Norte   12       5           340            100            1     2
## 2  760 Zona Norte    2       5           380            100            1     3
## 3  789 Zona Norte    1       5           300            100            1     2
## 4  946 Zona Norte   12       5           325            100            1     3
## 5 1093 Zona Norte    2       5           240            100            1     2
##   Habitaciones        Tipo           Barrio cordenada_longitud
## 1            3 Apartamento         la flora          -76.49995
## 2            2 Apartamento       zona norte          -76.50274
## 3            3 Apartamento         la flora          -76.50300
## 4            3 Apartamento         la flora          -76.50623
## 5            3 Apartamento prados del norte          -76.51000
##   Cordenada_latitud
## 1           3.48347
## 2           3.42428
## 3           3.48200
## 4           3.40303
## 5           3.48000

require(leaflet)
leaflet() %>% addCircleMarkers(lng = db_filtrada3$cordenada_longitud,lat = db_filtrada3$Cordenada_latitud,radius = 0.3,color = "black") %>% 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.

library(readxl)
datos = read_excel("C:/Users/User/Desktop/data arboles.xlsx")

datos$peso = as.numeric(datos$peso)
datos$diametro = as.numeric(datos$diametro)
datos$altura = as.numeric(datos$altura)

head(datos)

## # 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.9

plot(datos$altura, datos$peso, 
     pch=16, 
     las=1,
     ylab = "peso del árbol",
     xlab = "altura del árbol")
grid()

cor(datos$altura, datos$peso)

## [1] 0.8582009

modelo1=lm(peso ~ altura, data= datos)
summary(modelo1)

## 
## Call:
## lm(formula = peso ~ altura, data = datos)
## 
## 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-16

Se identifica que hay una relación fuerte entre la variable altura y peso, es decir, que el peso del árbol depende en un 73% de la altura.

cor(datos[3:5])

##               peso diametro    altura
## peso     1.0000000 0.908123 0.8582009
## diametro 0.9081230 1.000000 0.9355360
## altura   0.8582009 0.935536 1.0000000

plot(datos[3:5])

modelo2 = lm( peso ~ altura + diametro, data = datos)
summary(modelo2)

## 
## Call:
## lm(formula = peso ~ altura + diametro, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3083 -2.5121  0.1608  2.0088 11.7446 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.1205     1.4305  -6.376 8.44e-09 ***
## altura        0.3132     0.5751   0.544    0.587    
## diametro      4.7395     0.7128   6.649 2.49e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.449 on 87 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8213 
## F-statistic: 205.5 on 2 and 87 DF,  p-value: < 2.2e-16

#Comparamos los dos modelos

library(memisc) 
memisc::mtable(modelo1, modelo2)

## 
## Calls:
## modelo1: lm(formula = peso ~ altura, data = datos)
## modelo2: lm(formula = peso ~ altura + diametro, data = datos)
## 
## =====================================
##                 modelo1    modelo2   
## -------------------------------------
##   (Intercept)  -7.046***  -9.121***  
##                (1.705)    (1.430)    
##   altura        3.891***   0.313     
##                (0.248)    (0.575)    
##   diametro                 4.739***  
##                           (0.713)    
## -------------------------------------
##   R-squared     0.737      0.825     
##   N            90         90         
## =====================================
##   Significance: *** = p < 0.001;   
##                 ** = p < 0.01;   
##                 * = p < 0.05

En este comparativo se evidencia que el Modelo 2 es mas representativo con un 82% y este es el seleccionado.

Ahora definimos variables:

modelo3= step(modelo2)

## Start:  AIC=225.79
## peso ~ altura + diametro
## 
##            Df Sum of Sq    RSS    AIC
## - altura    1      3.53 1038.2 224.09
## <none>                  1034.7 225.79
## - diametro  1    525.74 1560.5 260.76
## 
## Step:  AIC=224.09
## peso ~ diametro
## 
##            Df Sum of Sq    RSS    AIC
## <none>                  1038.2 224.09
## - diametro  1      4884 5922.2 378.80

summary(modelo3)

## 
## Call:
## lm(formula = peso ~ diametro, data = datos)
## 
## 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-16

Siempre buscando la mejor reprsentación del modelo a la realidad, buscamos alternativas y una de ellas es hacer una transformación Log

modelo4= lm(log(peso) ~ diametro, data = datos)
summary(modelo4)

## 
## Call:
## lm(formula = log(peso) ~ diametro, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27395 -0.10180 -0.00328  0.10073  0.33742 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.32798    0.05977   22.22   <2e-16 ***
## diametro     0.27818    0.01061   26.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1453 on 88 degrees of freedom
## Multiple R-squared:  0.8865, Adjusted R-squared:  0.8852 
## F-statistic: 687.6 on 1 and 88 DF,  p-value: < 2.2e-16

#Comparamos modelos

memisc::mtable(modelo3, modelo4)

## 
## Calls:
## modelo3: lm(formula = peso ~ diametro, data = datos)
## modelo4: lm(formula = log(peso) ~ diametro, data = datos)
## 
## =====================================
##                 modelo3    modelo4   
##                ---------  ---------  
##                  peso     log(peso)  
## -------------------------------------
##   (Intercept)  -9.020***   1.328***  
##                (1.413)    (0.060)    
##   diametro      5.103***   0.278***  
##                (0.251)    (0.011)    
## -------------------------------------
##   R-squared     0.825      0.887     
##   N            90         90         
## =====================================
##   Significance: *** = p < 0.001;   
##                 ** = p < 0.01;   
##                 * = p < 0.05

Ahora se observa que al hacer la transformación, aumenta la representación del modelo a la realidad a un 89%, es por eso, que este es el modelo ideal para trabajar.

VALIDACIÓN CRUZADA:

Se realizará un escenario 80-20(modelar-validar)

#Paso 1, Segmentar los datos

id_modelar=sample(1:200,size=160)
datos_modelar=datos[id_modelar,]
datos_validar=datos[-id_modelar,]


#Paso 2: Estimar el modelo Set Modelar

mod_datos_modelar=lm(peso ~ diametro, data=datos_modelar)

#Paso 3: Predecir set de validación

peso_pred=predict(mod_datos_modelar,list(diametro=datos_validar$diametro))

#Paso 4: Comparar modelo con realidad

peso_real = datos_validar$diametro
error= peso_real-peso_pred

#Paso 5: - Calcular Indicador de evaluación de la predicción
MAE=mean(abs(error))
MAE

## [1] 13.1112

La validación cruzada en un primer paso segmentamos los datos dejando 80% para el modelo y 20% aleatorios para validar, luego se ajusta el modelo con el 80%. Finalmente predecirmos el peso del 20% de arboles y se comparan los resultados contra los reales por medio de la metrica MAE, el cual nos da alrededor de 13, es decir, este modelo tiene un error de predicción alrededor de 13 toneladas en este caso de los arboles.

Informe de aplicación de regresión lineal

Santiago Garcia Cifuentes

17/8/2022

https://rpubs.com/sgcifuentes/935217