Taller 1 Regresión Lineal
Julieth Cerón
2022-08-15
INFORME NO.1 DE APLICACIÓN REGRESIÓN LINEAL SIMPLE
PUNTO 1. PREDICCIÓN DE LAS ACCIONES
Analizar el comportamiento de los precios de las Acciones de Ecopetrol (p_ecop) según la variación del precio del barril de petróleo WTI (p_wti)producido en Colombia. Los datos son los siguientes:
fech= 1:18
p_ecop=c(1090, 1170, 1160, 1230, 1155, 1165, 1205, 1170, 1150, 1130, 1110, 1105, 1085, 1060, 1035, 1015, 955, 961)
p_wti=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)
data1= data.frame(fech, p_ecop, p_wti)
attach(data1)## The following objects are masked _by_ .GlobalEnv:
##
## fech, p_ecop, p_wtihead(data1, n = 18)## fech p_ecop p_wti
## 1 1 1090 35.62
## 2 2 1170 36.31
## 3 3 1160 37.35
## 4 4 1230 34.95
## 5 5 1155 34.53
## 6 6 1165 35.81
## 7 7 1205 36.14
## 8 8 1170 37.50
## 9 9 1150 37.80
## 10 10 1130 36.81
## 11 11 1110 37.87
## 12 12 1105 37.04
## 13 13 1085 36.76
## 14 14 1060 35.97
## 15 15 1035 33.97
## 16 16 1015 33.27
## 17 17 955 31.41
## 18 18 961 30.44A continuación se presenta de forma descriptiva el comportamiento de las variables de estudio:
Indicadores de variabilidad del precio de las acciones de ecopetrol
En la tabla siguiente se observa que en promedio las acciones de Ecopetrol cuestan $1.108 con una desviación estandar de $78.42 y una variación de $6.150.
library(moments)
Indicadores_pecop=data.frame("Media"=mean(p_ecop),"Mediana"=median(p_ecop),"Desvest"=sd(p_ecop),"Varianza"=var(p_ecop))
Indicadores_pecop## Media Mediana Desvest Varianza
## 1 1108.389 1120 78.42354 6150.252Indicadores de variabilidad del precio del barril de petróleo
En la tabla siguiente se observa que valor medio de un barril de petróleo esta en 35.53 dólares con una desviación estandar de 2.12 dólares y una varianza de 4.48 dólares.
Indicadores_pwti=data.frame("Media"=mean(p_wti),"Mediana"=median(p_wti),"Desvest"=sd(p_wti),"Varianza"=var(p_wti))
Indicadores_pwti## Media Mediana Desvest Varianza
## 1 35.53056 36.055 2.118183 4.4867Para observar la relación de estas dos variables de análisis, se grafica a continuación un diagrama de dispersión y se cálcula el coeficiente de correlación de Pearson para ver si es posible explicar la variabilidad del precio de las acciones de Ecopetrol (p_ecop) mediante la variación del precio del barril de petróleo (p_wti):
library(ggplot2)
ggplot(data1, aes(x = p_wti, y = p_ecop)) + geom_point(colour="purple") + geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
cor(p_ecop,p_wti)## [1] 0.7074373En el gráfico anterior se observa una relación entre esta dos variables, que inicialmente pareciera no ser necesariamente lineal. Sin embargo, se encuentra que el coeficiente de correlación de Pearson es igual al 70.74%, con lo que se podría considerar que inicialmente si se podría generar la regresión lineal entre ambas variables.
a.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.
Modelo lineal propuesto:Precio de las acciones de Ecopetrol= β0 + β1 Precio del petróleo WTI + e
lm=lm(p_ecop~p_wti)
resumen=summary(lm)
resumen##
## Call:
## lm(formula = p_ecop ~ p_wti)
##
## 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
## p_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.001024La estimación de la recta (ecuación de regresión) quedaría de la siguiente forma:
Precio de las acciones de Ecopetrol= 177.77 + 26.19 * Precio del petróleo WTI + e
El coeficiente de determinación (R2) representa la proporción de la variabilidad de Y que es posible explicar a travez de x. En este caso el modelo construido explica sólo el 50.05% de las variaciones del precio de las acciones de Ecopetrol a través del precio del petróleo WTI.
**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**
Para analizar la significacia del modelo, se debe validar que el β^1 sea significativo, por tanto, se evaluan las siguientes hipótesis:
H0:β1=0 H1:β1≠0
resumen$coefficients[2, ]## Estimate Std. Error t value Pr(>|t|)
## 26.192134618 6.541913749 4.003741966 0.001023938Al realizar la prueba de significancia observamos que el valor-p asociado al coeficiente del precio del barril de petróleo (p_wti) es de 0.001, lo que nos dice que contrastando este resultado con un nivel de significacia del 0.05, no se cuenta con evidencia estadistica suficiente para suponer que el β1 estimado no sea significativo o en otras palabras, el valor de la pendiente es un predictor significativo para la ecuación del precio de la acción de Ecopetrol.
c. Interprete los coeficientes del modelo propuesto en “a)”
De acuerdo al punto a, la interpretación de los betas sería:
β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.
β1: Este nos quiere decir que por cada dolar adicional que incremente o sume el precio del barril de petróleo el precio de las acciones de Ecopetrol incrementará en $26.192.
d. Haga un análisis de los residuos. ¿Qué supuesto no se cumple?
Supuesto de Linealidad: Si se cumple ya que el valor esperado de los residuos es igual a cero.
library(ggplot2)
yest = predict(lm)
ggplot(data1,aes(x = p_wti, y = yest)) +
geom_point() +
geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm$residuals),3)))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
mean(lm$residuals)## [1] -5.527407e-15Supuesto de Normalidad:
Como p-value = 0.04276 < α = 0.05, se rechaza la hipótesis nula, quiere decir que, los errores no se distibuyen de forma normal, el supuesto no se cumple.
Las hipotesis a contrastar son:
H0: Los errores presentan una distribución normal
H1: Los errores no presentan una distribución normal
shapiro.test(lm$residuals)##
## Shapiro-Wilk normality test
##
## data: lm$residuals
## W = 0.89259, p-value = 0.04276En la gráfica siguiente de residuales vs ajustados se evidencia que los residuos no se comportan de forma aleatoria como se esperaría. En conclusión no se cumplen los supuestos relacionados con la normalidad y la aleatoriedad de los residuos.
errorlm=lm$residuals
par(mfrow=c(2,2))
plot(lm)
Supuesto de Homocedasticidad de Varianza (Breush Pagan)
Como $ P-value $ es mayor a 0.05 (nivel de significancia escogido), no se rechaza H0, entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
H0:Los residuales se distribuyen con la misma varianza
Ha:Los residuales NO se distribuyen con la misma varianza
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(lm)##
## studentized Breusch-Pagan test
##
## data: lm
## BP = 0.029563, df = 1, p-value = 0.8635Supuesto de Independencia o autocorrelación
Como $ P-value $ (0.001) es menor a 0.05 (nivel de significancia escogido),se rechaza H0, entonces se podría pensar que los errores podrian estar autocorrelacionados.Es decir que no se cumple con el supuesto de independencia de los errores.
Las hipotesis planteadas son:
H0: Los errores no presentan autocorrelación
H1: Los errores presentan autocorrelación
library(lmtest)
lmtest::dwtest(lm)##
## Durbin-Watson test
##
## data: lm
## DW = 0.74504, p-value = 0.0004666
## alternative hypothesis: true autocorrelation is greater than 0library(lmtest)
dw = dwtest(lm, alternative = "two.sided", iterations = 1000)
plot(lm$residuals, main = paste("dw P-Value = ",round(dw$p.value,3)))
e. Concluya sobre la validez del modelo propuesto en a)
Analizando el modelo propuesto, se observa que al utilizar unicamente el precio del petroleo no se puede llegar a predecir correctamente el precio de las acciones de ecopetrol dado que el modelo propuesto solo explica el 50% de la informacion el precio de las acciones, además el modelo no cumple con el supuesto de la normalidad, aletoriedad e independencia de los errores, por tanto, se recomienda explorar una transformación del modelo para validar si hay un mejor comportamiento del mismo.
PUNTO 2.Los siguientes datos corresponden a la INFLACION (inf) y al SALARIO MINIMO MENSUAL LEGAL VIGENTE (smml) desde el año 1999 para Colombia.
fec = 1999:2015
inf=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)
smml=c(236460, 260100, 286000, 309000, 332000, 358000, 381500, 408000, 433700, 461500, 496900, 515000, 535600, 566700, 589500, 616027, 644350)
data2 = data.frame(fec, inf, smml)
attach(data2)## The following objects are masked _by_ .GlobalEnv:
##
## fec, inf, smmlhead(data2, n = 17)## fec inf smml
## 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 358000
## 7 2005 4.85 381500
## 8 2006 4.48 408000
## 9 2007 5.69 433700
## 10 2008 7.67 461500
## 11 2009 2.00 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 644350La 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 (smml) y como variable independiente INFLACION (inf, obtenga un modelo de regresión lineal simple y resuelva:
A continuación se presenta de forma descriptiva el comportamiento de las variables de estudio:
Indicadores de variabilidad de la inflación (inf)
En la tabla siguiente se observa que en promedio la inflación en Colombia fue de 5.35% entre 1999 y 2015 con una desviación estandar de 2.32% y una varianza de 5,37%.
library(moments)
Indicadores_inf=data.frame("Media"=mean(inf),"Mediana"=median(inf),"Desvest"=sd(inf),"Varianza"=var(inf))
Indicadores_inf## Media Mediana Desvest Varianza
## 1 5.353529 5.5 2.318253 5.374299Indicadores de variabilidad del salario mínimo (smml)
En la tabla siguiente se observa queel valor promedio del salario mínimo en Colombia entre 1999 y 2015 fue de $ 437.078.7 con una desviación estandar de $ 129.183.
library(moments)
Indicadores_smml=data.frame("Media"=mean(smml),"Mediana"=median(smml),"Desvest"=sd(smml),"Varianza"=var(smml))
Indicadores_smml## Media Mediana Desvest Varianza
## 1 437078.6 433700 129182.6 16688136605Para observar la relación de estas dos variables de análisis, se grafica a continuación un diagrama de dispersión y se cálcula el coeficiente de correlación de Pearson para ver si es posible explicar la variabilidad del salario mínimo (smml) mediante la variación de la inflación (inf):
En el gráfico siguiente 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.Al calcular el coeficiente de correlación de Pearson se obtiene un valor de -70.87%, es decir que entre ambas varibles efectivamente existe una relación lineal negativa débil.
library(ggplot2)
ggplot(data2, aes(x = inf, y = smml)) + geom_point(colour="purple") + geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
cor(smml,inf)## [1] -0.7086581a. Escriba la ecuación del modelo de regresión lineal simple
Salariomínimo=β0+β1Inflación+ϵ
lm2 = lm(smml~inf)
resumen2=summary(lm2)
resumen2##
## Call:
## lm(formula = smml ~ inf)
##
## 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 ***
## inf -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.00145La ecuación de la recta quedaría expresada de la siguiente forma:
Salariomínimo=648.486 − 39.489∗Inflacion + e
Al analizar la significacia del modelo, lo hacemos validando que el β^1 sea significativo, por tanto la hipotesis a evaluar son las siguientes:
H0:β1=0 H1:β1≠0
Al realizar la prueba de significancia observamos que el valor-p asociado al coeficiente de la variable inflación es de 0.00145, lo que nos dice que contrastando este resultado con un nivel de significacia del 0.05, no se cuenta con evidencia estadistica suficiente para suponer que el β1 estimado no sea significativo.
b.Plantee y valide las hipótesis correspondientes a la linealidad general del modelo propuesto en a)
Para analizar la significacia del modelo, se debe validar que el β^1 sea significativo, por tanto, se evaluan las siguientes hipótesis:
H0:β1=0 H1:β1≠0
Al realizar la prueba de significancia observamos que el valor-p asociado al coeficiente del precio del barril de petróleo (p_wti) es de 0.00145, lo que nos dice que contrastando este resultado con un nivel de significacia del 0.05, no se cuenta con evidencia estadistica suficiente para suponer que el β1 estimado no sea significativo o en otras palabras, el valor de la pendiente es un predictor significativo para la ecuación del salario mínimo en Colombia.
resumen2$coefficients[2, ]## Estimate Std. Error t value Pr(>|t|)
## -3.948933e+04 1.015141e+04 -3.890033e+00 1.450420e-03Supuesto de Linealidad:
H0: ß0 igual a 0 H1: ß0 diferente de 0
Se cumple el supuesto de linealidad ya que el valor promedio d elos residuos es cero:
mean(lm2$residuals)## [1] -1.491304e-12t.test(lm2$residuals, mu=0)##
## One Sample t-test
##
## data: lm2$residuals
## 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-12p-value = 1 > α = 0.05, se tiene evidencia suficiente para no rechazar que H0: ß0 = 0, por tanto el supuesto se cumple.
yest2 = predict(lm2)
ggplot(data2,aes(x = inf, y = yest2)) + geom_point() + geom_smooth() +
labs(title = paste("Media residuos =", round(mean(lm2$residuals),3)))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
c. Indique e interprete el coeficiente de correlación del modelo propuesto en a)
Al calcular el coeficiente de correlación de Pearson se obtiene un valor de -70.87%, es decir que entre el salario mínimo y la inflación existe una relación lineal negativa débil.
cor(smml,inf)## [1] -0.7086581d. Interprete cada uno de los coeficientes del modelo propuesto en a)
β0: En el caso dado donde la inflación llegase a ser 0, el valor del salario mínimo sería de $ 648.485,93.
β1: Por cada unidad porcentual adicional que incremente la inflación, el salario mínimo disminuirá en $39.489,33.
e. Construya una gráfica de residuales y haga un análisis cualitativo de los supuestos del modelo propuesto en a)
En las gráficas siguientes se evidencia que los errores no se distribuyen de forma aleatoria y que los datos no se ajustan perfectamente a la línea de normalidad. Ante ello, es recomendable transformar el modelo de regresion lineal para mejorarlo y corregir estos inconvenientes.
e2=lm2$residuals
par(mfrow=c(2,2))
plot(lm2)
Supuesto de Normalidad:
La prueba de Shapiro - Wilk arroja que el p_value es 0.001407 < 0.05, por lo tanto, se rechaza la hipotesis nula y se concluye que los residuos no siguen una distribución normal.
shapiro.test(lm2$residuals)##
## Shapiro-Wilk normality test
##
## data: lm2$residuals
## W = 0.78826, p-value = 0.001407Supuesto de Varianza constante:
Según la prueba Breusch-Pagan, el p-value 0.433 > 0.05, sugiere que no se rechaza Ho, por lo que se evidencia que la varianza de los residuos es homocedastica.
Ho = Homocedasticidad (los residuos se distibuyen con la misma varianza) Ha = Heterocedasticidad (los residuos no se distribuyen con la misma varianza),
library(lmtest)
bptest(lm2)##
## studentized Breusch-Pagan test
##
## data: lm2
## BP = 0.61478, df = 1, p-value = 0.433Supuesto de independencia de los errores
Como $ P-value $ (0.0002) es menor a 0.05 (nivel de significancia escogido),se rechaza H0, entonces se podría pensar que los errores podrian estar autocorrelacionados.Es decir que no se cumple con el supuesto de independencia de los errores.
Las hipotesis planteadas son:
H0: Los errores no presentan autocorrelación
H1: Los errores presentan autocorrelación
library(lmtest)
lmtest::dwtest(lm2)##
## Durbin-Watson test
##
## data: lm2
## DW = 0.68432, p-value = 0.0002714
## alternative hypothesis: true autocorrelation is greater than 0f.Comente sobre la conveniencia de usar el modelo propuesto en a) para predecir el SMLM para Colombia.
El modelo propuesto en el punto a) no seria conveniente para predecir el valor del SMML ya que registra un valor de R2 muy bajo. Además, invalida el supuesto de normalidad, aletoriedad e independencia de los errores. Se debe destacar que,el salario mínimo puede estar influenciado por otros factores como, la Producción Total de los Factores, acuerdos bilaterales entre Gobierno y Central de trabajadores, entre otros y no unicamente de la inflacion. Adicionalmente, la inconsistencia sobre el comportamiento normal de las variables llama la atención dado que se esperaria que a mayor inflación el salario minimo fuera mayor. Una posible alternativa seria tomar como variable dependiente la variación % del salario mínimo y no su valor neto real.
PUNTO 3. MODELO DE REGRESIÓN MÚLTIPLE-PRECIOS DE LAS VIVIENDAS
Base de datos original y sin filtros:
library(readxl)
library(ggplot2)
datos_viviendas <- read_excel("D:/Maestria Ciencia de Datos/SEMESTRE II/Modelos Estadisticos/Datos_Vivienda.xlsx")
attach(datos_viviendas)
head(datos_viviendas,n = c(6,7))## # A tibble: 6 × 7
## 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 4a. 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?).
#Realizamos el filtro requerido
filtro = which ( datos_viviendas$Tipo == "Apartamento" & datos_viviendas$Zona == "Zona Norte" & datos_viviendas$precio_millon < 500 & datos_viviendas$Area_contruida < 300)
#Base de datos filtrada
viviendas_fil = datos_viviendas[filtro,]
# Los datos de la variable parqueadero se pasan a tipo numérico
viviendas_filt <- transform(viviendas_fil, parqueaderos = as.numeric(parqueaderos))## Warning in eval(substitute(list(...)), `_data`, parent.frame()): NAs
## introducidos por coerción#Se reemplazan los NA y queda la base final viviendas_filt
viviendas_filt[is.na(viviendas_filt)] <- 0
head(viviendas_filt, n = c(6,7))## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 1 Zona Norte 2 3 135 56 1 1
## 2 Zona Norte NA 3 78 54 2 1
## 3 Zona Norte NA 5 340 106 2 2
## 4 Zona Norte 1 3 135 103 1 3
## 5 Zona Norte 1 3 75 54 1 2
## 6 Zona Norte NA 4 175 77 1 2dim(datos_viviendas)[1]## [1] 8322dim(viviendas_filt)[1]## [1] 1077Inicialmente se contaba con 8322 registros, después de realizar los 4 filtros sugeridos, quedamos con 1077 registros.
En el mapa siguiente donde se referencia todos los apartamentos ubicados en la zona norte de la ciudad,ademas con precios inferiores a los 500 millones y areas menores a 300 m2, se puede observar que muchos de los puntos se encuentran referenciados en la zona sur, incluso algunos en la zona oeste y oriente, esto se debe a que en la base de datos original a muchos apartamentos que quedan en la zona sur, oeste y oriente le han reistrado coordenada longitud y latitud norte, en este caso como el modelo grafica las ubicaciones respecto a la longitud y latitud entonces se esta presentando el error.
require(leaflet)## Loading required package: leafletleaflet() %>% addCircleMarkers(lng= viviendas_filt$cordenada_longitud, lat=viviendas_filt$Cordenada_latitud, radius = 4, color = "purple") %>% addTiles()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.
Precio de las viviendas VS Área construida
En la figura siguiente se observa que si existe una relación positiva entre el precio y el area construida, es decir que, mientras el area construida sea mayor el precio será más alto:
Fig1=ggplot(viviendas_filt,aes(y=precio_millon,x=Area_contruida)) + geom_point(colour = "purple", size = 2) + geom_smooth()
Fig1## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
library(plotly)##
## Attaching package: 'plotly'## 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':
##
## layoutgraf1 <- plot_ly(data = viviendas_filt, x =~viviendas_filt$Area_contruida , y = ~viviendas_filt$precio_millon ,
marker = list(size = 10,
color = 'purple',
line = list(color = 'black',
width = 2)))
graf1 <- graf1 %>% layout(title = 'Precio Vivienda vs Area Construida',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
graf1## 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-modePrecio de las viviendas VS Estrato
En el gráfico siguiente se observa una relación directa entre las dos variables, es decir, entre mayor sea el estrato, mayor es el precio de la vivienda.
require(ggplot2)
require(plotly)
Fig2= ggplot(viviendas_filt = viviendas_filt, mapping = aes(x= viviendas_filt$Estrato, y=viviendas_filt$precio_millon))+ geom_point(colour="purple")+geom_smooth()
ggplotly(Fig2)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.Precio de las viviendas VS Parqueaderos
El gráfico siguiente no da cuenta de una relación existente entre el número de parqueaderos vs el precio de la vivienda.
require(ggplot2)
require(plotly)
Fig3= ggplot(viviendas_filt = viviendas_filt, mapping = aes(x= viviendas_filt$parqueaderos, y=viviendas_filt$precio_millon))+ geom_point(colour="purple")+geom_smooth()
ggplotly(Fig3)## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Computation failed in `stat_smooth()`:
## x has insufficient unique values to support 10 knots: reduce k.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).
La ecuación de regresión corresponde a:
Preciomillon = −157.38220 + 0.94938 (Area_contruida) + 68.99436 (Estrato) + 22.64906 (parqueaderos)
El valor del R2 = 0.7623 significa que las variables: area construida, estrato y parquaderos explican el 76% de la variabilidad del precio del apartamento.
Significancia ßo,ß1,ß2 y ß3 son estadisticamente significativos a un nivel de significancia de 0.05.
Hipótesis para ßo: Ho:βo=0,Ha:β0≠0
p−value=2e−16∗∗∗
Ya que p-value < 0.05, se rechaza H0: ß0 = 0 Hipótesis para ß1: H0:β1=0,Ha:β1≠0
p−value=2e−16∗∗∗
Ya que p-value < 0.05, se rechaza Ho: ß2 = 0 Hipótesis para ß2: Ho:β2=0,Ha:β2≠0
p−value=2e−16∗∗∗
ya que p-value < 0.05, se rechaza Ho: ß2 = 0 Hipótesis para ß3: Ho:β3=0,Ha:β3≠0
p−value=3.24e−16∗∗∗
Ya que p-value < 0.05, se rechaza Ho: ß3 = 0
Interpretación coeficientes:
ß1: Si el área construida de la vivienda aumenta en 1mts2, se espera que el precio de la vivienda aumente en 0,94 millones. ß2:A mayor estrato, el precio de la vivienda aumente en 22.64 millones. * ß3:Por cada cada parqueadero adicional, se espera que el precio la vivienda aumente en 15.55 millones *El coeficiente R2 ajustado de 0.76 indica que el precio de la vivienda es explicado por el modelo en un 76%
modelo_viviendas = lm(formula= precio_millon ~ Area_contruida + Estrato + parqueaderos, data = viviendas_filt)
summary(modelo_viviendas)##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Estrato + parqueaderos,
## data = viviendas_filt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -216.571 -31.564 -1.213 27.889 224.053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -157.38220 7.71190 -20.408 < 2e-16 ***
## Area_contruida 0.94938 0.06054 15.682 < 2e-16 ***
## Estrato 68.99436 2.26623 30.445 < 2e-16 ***
## parqueaderos 22.64906 2.73064 8.294 3.24e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.75 on 1073 degrees of freedom
## Multiple R-squared: 0.7629, Adjusted R-squared: 0.7623
## F-statistic: 1151 on 3 and 1073 DF, p-value: < 2.2e-16d. 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).
En el gráfico de residuales V.S valores ajustados se observa que hay 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.
par(mfrow = c(2, 2))
plot(modelo_viviendas)
Supuesto de normalidad
Como p-value = 3.509e-08 < α = 0.05, se rechaza la hipótesis nula. COn esto se concluye que los errores no se distibuyen de forma normal, el supuesto no se cumple. Ademásla gráfica de Normal QQ también valida este comportamiento.
shapiro.test(modelo_viviendas$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_viviendas$residuals
## W = 0.98702, p-value = 3.509e-08Supuesto de Linealidad: Si se cumple ya que el valor esperado de los residuos es igual a cero.
summary(modelo_viviendas$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -216.571 -31.564 -1.213 0.000 27.889 224.053t.test(modelo_viviendas$residuals, mu=0)##
## One Sample t-test
##
## data: modelo_viviendas$residuals
## t = 6.246e-16, df = 1076, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -3.209419 3.209419
## sample estimates:
## mean of x
## 1.021631e-15p-value = 1 > α = 0.05, se tiene evidencia suficiente para no rechazar que Ho: ßo = 0, por tanto el supuesto se cumple.
Supuesto de varianza constante
H0:Los residuales se distribuyen con la misma varianza
Ha:Los residuales NO se distribuyen con la misma varianza
library(lmtest)
bptest(modelo_viviendas)##
## studentized Breusch-Pagan test
##
## data: modelo_viviendas
## BP = 156.86, df = 3, p-value < 2.2e-16el valor p de la prueba es menor que el nivel de significancia (es decir, α = .05), entonces rechazamos la hipótesis nula y concluimos que la heterocedasticidad está presente en el modelo de regresión.
Supuesto de autocorrelación
H0:No existe correlación entre los errores
Ha:Existe correlación entre los errores
Como p-value = 2.992e-05 < α = 0.05, se rechaza la hipótesis nula, quiere decir que, los errores no son indepentientes, el supuesto no se cumple.
lmtest::dwtest(modelo_viviendas)##
## Durbin-Watson test
##
## data: modelo_viviendas
## DW = 1.757, p-value = 2.992e-05
## alternative hypothesis: true autocorrelation is greater than 0En este caso se viola el supuesto de normalidad, homocedasticidad, autocorrelación, por ende se podría pensar en realizar alguna transformación sobre la variable respuesta o sobre las variables explicativas (trabajar en logaritmos).
e. Con el modelo identificado predecir el precio de un apartamento con 100 mt2, de estrato 4 y con parqueadero. ¿Si este apartamento lo están ofreciendo en 450 millones cual seria su opinión con base en el resultado del modelo considera que es una buena oferta?
modelo_viviendas$coefficients## (Intercept) Area_contruida Estrato parqueaderos
## -157.3822027 0.9493794 68.9943565 22.6490604De esta forma, la ecuación nos queda de la siguiente manera:
precio_pred=-157.3822027+0.9493794*100+ 68.994*4 + 22.6490604*1
precio_pred## [1] 236.1808#con un parqueadero:
predict(modelo_viviendas,list(Area_contruida=100, Estrato = 4, parqueaderos = 1), interval = "confidence" )## fit lwr upr
## 1 236.1822 232.2449 240.1195Considerando los valores obtenidos del modelo, no es buena opción comprar el apartamento que estan ofreciendo por 450 millones, dado que el modelo registra un valor de compra aproximado de 236 millones.
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 encuenta 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
Aptos = which ( viviendas_filt$Area_contruida >100 & viviendas_filt$Estrato == 4 & viviendas_filt$parqueaderos >=1 & viviendas_filt$precio_millon <= 400 & viviendas_filt$precio_millon >= 236.1822)
#datos filtrados:
viviendas_filt2 = viviendas_filt[Aptos,]
min(viviendas_filt2$precio_millon)## [1] 240max(viviendas_filt2$Area_contruida)## [1] 287Considerando que la persona tiene un crédito preaprobado por un valor de 400 millones se sugiere al comprador ofertas de apartamento en zona norte, estrato 4 y con posibilidad de 1 o 2 parqueaderos y área construida mayor a 100mt2. A continuación se relacionan as ofertas posibles de adquirir en el barrio la flora o versalles, zonas de la ciudad de Cali que se encuentran muy bien valoradas por compradores:
Aptos1 = which((viviendas_filt2$precio_millon == 240 | viviendas_filt2$Area_contruida == 287 | viviendas_filt2$Habitaciones ==4 | viviendas_filt2$Banos ==3) & (viviendas_filt2$Barrio =="la flora" | viviendas_filt2$Barrio == "versalles") )
viviendas_Filt5 = viviendas_filt2[Aptos1,]
head(viviendas_Filt5,n = c(7,7))## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## 418 Zona Norte 4 4 380 123 1 3
## 755 Zona Norte NA 4 265 125 2 3
## 788 Zona Norte 2 4 380 126 2 3
## 800 Zona Norte 6 4 270 152 1 3
## 844 Zona Norte 2 4 240 103 1 2
## 893 Zona Norte NA 4 270 111 1 3
## 921 Zona Norte 7 4 300 126 2 4require(leaflet)
leaflet() %>% addCircleMarkers(lng= viviendas_Filt5$cordenada_longitud, lat=viviendas_Filt5$Cordenada_latitud,radius = 8, color = "purple") %>% addTiles()PUNTO 4. MODELO DE REGRESIÓN MÚLTIPLE
a. 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_arboles <- read_excel("D:/Maestria Ciencia de Datos/SEMESTRE II/Modelos Estadisticos/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)
head(datos_arboles,n = c(6,7))## # A tibble: 6 × 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.9Las variables que contiene la base de datos son: la finca y el tipo de genotipo como variables cualitativas y el peso, diametro, y altura de los arboles como variables cuantitativas . A continuaciòn se hara la construcciòn de dos modelos de regresiòn simple usando las variables continuas y al final se logrará determinar 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 0ANÁLISIS BIVARIADO
Peso VS Altura
En el gráfico siguiente se observa que hay una relación lineal positiva entre la altura y el peso del árbol. A mayor altura, mayor peso del árbol.
library(plotly)
Fig4=ggplot(datos_arboles,aes(y=peso,x=altura)) + geom_point(colour = "purple", size = 2) + geom_smooth()
Fig4## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Peso VS Diámetro
En el gráfico siguiente se observa que hay una relación lineal positiva entre la altura y el diámetro del árbol. A mayor diámetro, mayor peso del árbol.
library(plotly)
Fig5=ggplot(datos_arboles,aes(y=peso,x=diametro)) + geom_point(colour = "purple", size = 2) + geom_smooth()
Fig5## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Modelo de regresión múltiple
modelo_arboles = lm(peso ~ diametro + altura, data = datos_arboles)
summary(modelo_arboles)##
## 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-16La ecuación de regresión corresponde a: Peso = -9.1205 + 4.7395(diametro) + 0.3132(altura)
Significancia del modelo
Hipótesis para ßo: Ho:βo=0,Ha:β0≠0 p−value=8.44e-09∗∗∗ Ya que p-value < 0.05, se rechaza Ho: ß0 = , es estadisticamente significativo.
Hipótesis para ß1: Ho:β1=0,Ha:β1≠0 p−value=2.49e-09∗∗∗ Ya que p-value < 0.05, se rechaza Ho: ß1 = ; es estadisticamente significativo.
Hipótesis para ß2: Ho:β2=0,Ha:β2≠0 p−value=0.587∗∗∗ ya que p-value > 0.05, no se rechaza Ho: ß2 = 0;No es estadisticamente significativo.
Interpretación:
R cuadrado es igual a 0.8253 significa que la variable: diametro y altura explican en conjunto el 82,5% el peso del árbol.
Se pueden interpretar los βetas de la siguiente forma:
*B1: Si el diámetro del árbol aumenta en una unidad, se espera que el peso del árbol aumente en 4.73 unidades.
*B2( en caso de que fuera significativo):Por cada metro adicional en la altura del árbol, se esperaría que el peso del árbol aumente en 0.31 unidades.
Validación supuestos
par(mfrow=c(2,2))
plot(modelo_arboles)
e_medio=mean(residuals(modelo_arboles))
e_medio## [1] 1.032498e-16El valor esperado (promedio) de los residuales es igual a 0.
Supuesto de normalidad
norm_arb=shapiro.test(modelo_arboles$residuals)
norm_arb##
## Shapiro-Wilk normality test
##
## data: modelo_arboles$residuals
## W = 0.95745, p-value = 0.004966Com P−value es menor a a 0.05 (nivel de significancia escogido), se rechaza H0, entonces podría pensar que los errores no siguen una distribución normal.
Supuesto de Homocedasticidad de Varianza (Breush Pagan)
H0:Homocedasticidad del error
Ha:Heterocedasticidad del error
homocedasticidad_arb=bptest(modelo_arboles)
homocedasticidad_arb##
## studentized Breusch-Pagan test
##
## data: modelo_arboles
## BP = 14.32, df = 2, p-value = 0.0007772Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), no se rechaza H0, entonces se podría pensar que los errores cumplen con el supuesto de homocedasticidad.
Supuesto de Autocorrelación de los errores (Durbin-Watson)
H0: errores no autocorrelacionados
Ha:errores autocorrelacionados
lmtest::dwtest(modelo_arboles)##
## Durbin-Watson test
##
## data: modelo_arboles
## DW = 1.0481, p-value = 4.105e-07
## alternative hypothesis: true autocorrelation is greater than 0Como $ P-value $ es menor a 0.05 (nivel de significancia escogido), se rechaza H0, entonces se podría pensar que los errores podrian estar autocorrelacionados.
Validación cruzada:
Se utiliza una muestra del 80% de los datos para el modelo y el 20% para realizar las validaciones:
# Conjunto de 100 datos (80% para modelar y 20% para validar):
id_modelar= sample(1:90, size=72)
id_modelar## [1] 62 39 6 77 51 11 15 36 12 9 45 5 40 23 47 52 67 4 54 78 22 64 10 73 18
## [26] 31 56 82 21 65 37 41 16 8 89 87 44 17 50 74 59 57 72 26 7 48 2 68 79 32
## [51] 30 83 25 34 58 81 38 28 84 20 75 70 90 63 66 55 24 14 86 53 42 61# conjunto de datos con 72 datos para modelar:
arboles_modelar= datos_arboles[id_modelar,]
head(arboles_modelar)## # A tibble: 6 × 5
## finca mg peso diametro altura
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 FINCA_3 GENOTIPO_2 24.5 6 7.2
## 2 FINCA_2 GENOTIPO_1 20.8 6.2 6.8
## 3 FINCA_1 GENOTIPO_2 19.3 6.3 7.9
## 4 FINCA_3 GENOTIPO_1 13.0 4.2 4.8
## 5 FINCA_2 GENOTIPO_2 30.5 7.2 9.4
## 6 FINCA_1 GENOTIPO_1 7.03 2.5 3.7# Conjunto de datos con 18 datos para validar:
arboles_validar= datos_arboles[-id_modelar,]
head(arboles_validar)## # A tibble: 6 × 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 15.9 4.8 5.8
## 3 FINCA_1 GENOTIPO_1 7.47 2.2 3.5
## 4 FINCA_1 GENOTIPO_2 12 4.9 7
## 5 FINCA_1 GENOTIPO_2 9.89 4.3 6.3
## 6 FINCA_1 GENOTIPO_2 14.4 5.5 7.4# Estimar el Modelo del conjuto de datos de arboles_modelar
modelo_arboles1 = lm(formula= peso ~ diametro + altura, data = arboles_modelar)
# Predeccir el conjunto de datos de arboles_validar
predecir_peso = predict(modelo_arboles1, list(diametro=arboles_validar$diametro, altura=arboles_validar$altura))
# comparar pesos arrojados por modelo y los datos reales
peso_real = arboles_validar$peso
error = peso_real-predecir_peso
Tabla_Comparativa = data.frame(peso_real, predecir_peso, error)
head(Tabla_Comparativa, 5)## peso_real predecir_peso error
## 1 13.73 14.61001 -0.8800146
## 2 15.88 15.40913 0.4708725
## 3 7.47 2.21315 5.2568505
## 4 12.00 16.37215 -4.3721529
## 5 9.89 13.25758 -3.3675796#Mean Error Absoluto
MAE = mean(abs(error))
MAE## [1] 2.801006MAE: 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. MAE indica que la diferencia entre el peso real vs el peso que se predijo en promedio, es de 2.60 unidades de peso del arbol.
#RMSE
library(Metrics)
rmse(Tabla_Comparativa$peso_real, Tabla_Comparativa$predecir_peso) ## [1] 3.378745RMSE: 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. En otras palabras, la métrica RMSE indica que son los errores en unidades de masa en el peso del árbol, pero esta ultima medida tiene mayor impacto en los residuales más grandes, cuando son elevados al cuadrado. El valor para RMSE es pequeño, por lo que se puede decir que Los valores más bajos de RMSE indican un mejor ajuste.
Teniendo en cuenta las validaciones se puede indicar que el modelo se equivoca en un PORCENTAJE PEQUEÑO del peso promedio de acuerdo al MAE y la desviación estándar de la varianza inexplicada tambien es pequeña de acuerdo al RMSE.