Comenzamos leyendo la base de datos vivienda
## # A tibble: 8,322 × 12
## Zona piso Estrato precio_millon Area_contruida parqueaderos Banos
## <chr> <chr> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 Zona Sur 2 6 880 237 2 5
## 2 Zona Oeste 2 4 1200 800 3 6
## 3 Zona Sur 3 5 250 86 NA 2
## 4 Zona Sur NA 6 1280 346 4 6
## 5 Zona Sur 2 6 1300 600 4 7
## 6 Zona Sur 3 6 513 160 2 4
## 7 Zona Sur 2 6 870 490 3 6
## 8 Zona Sur 5 5 310 82.5 1 2
## 9 Zona Sur 9 4 240 80 1 2
## 10 Zona Sur 6 6 690 150 2 5
## # … with 8,312 more rows, and 5 more variables: Habitaciones <dbl>, Tipo <chr>,
## # Barrio <chr>, cordenada_longitud <dbl>, Cordenada_latitud <dbl>
Se procede a realizar la seleccion de las variables:
## # A tibble: 202 × 3
## precio_millon Area_contruida Barrio
## <dbl> <dbl> <chr>
## 1 1450 1200 el ingenio
## 2 290 100 el ingenio
## 3 360 99 el ingenio
## 4 550 197 el ingenio
## 5 410 136 el ingenio
## 6 390 198 el ingenio
## 7 300 147 el ingenio
## 8 250 97 el ingenio
## 9 370 104 el ingenio
## 10 1200 616 el ingenio
## # … with 192 more rows
se procede a realizar analisis descriptivo:
Observando el histograma de la variable precio millon, se puede decir que la mayoria de las viviendas ubicadas en el barrio el ingenio tienen un costo entre 0 y 500 millones de pesos; sin embargo tambien se observa que algunas viviendas pueden llegar a costar hasta mas de 1000 millones de pesos en el mismo barrio.
Observando el histograma de la variable area construida, se dice que la mayoria de las viviendas ubicadas en el barrio el ingenio, presentan un area entre 0 y los 200 metros cuadrados; tambien se pueden observar algunas viviendas que presentan hasta el doble o el triple del area en el mismo barrio.
Observando el grafico de dispersión, se dice que la mayoría de los datos de las viviendas en el barrio el ingenio, se encuentran concentrados en un rango de 0 a 200 millones de pesos junto con un area de 0 a 400 metros cuadrados.
Ahora se procede a realizar un analisis por medio de RLS:
##
## Call:
## lm(formula = precio_millon ~ Area_contruida, data = viviendanew)
##
## Residuals:
## Min 1Q Median 3Q Max
## -282.09 -61.15 -8.32 45.89 477.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 195.71496 11.63854 16.82 <2e-16 ***
## Area_contruida 1.21159 0.04632 26.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 105.2 on 200 degrees of freedom
## Multiple R-squared: 0.7738, Adjusted R-squared: 0.7727
## F-statistic: 684.2 on 1 and 200 DF, p-value: < 2.2e-16
Observando el modelo, se tiene que tanto el intercepto como la variable de area construida son significativos. Con respecto al coeficiente \(\beta_{1}\), se dice que por cada unidad de area construida va a aumentar en 1.21159 para la variable respuesta (precio millon). Ahora, observando el \(R^{2}\) del modelo indica que el modelo calculado explica el 77.38% de la variabilidad presente en la variable respuesta y el \(R^{2}\) ajustado fue de 77.27%.
Se procede a realizar la validación de los supuestos:
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 46.496, df = 1, p-value = 9.183e-12
Realizando el test de Breusch-pagan obtuvimos un valor p menor a 0.05 de tal manera que rechazamos la hipotesis nula y de esta manera no se cumple el supuesto de homocedasticidad.
###RESIDUOS
res.stud.base = studres(modelo)
modelo.fit.base= modelo$fitted.values
par(mfrow=c(1,2))
plot(modelo.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modelo.fit.base), col = 2)
plot(modelo.fit.base,abs(res.stud.base),
ylab='valor absoluto de los residuos estudentizados',
xlab='valores ajustados',main='(b)', col = "BLUE")
lines(lowess(abs(res.stud.base)~modelo.fit.base), col = 2)
Comprobando el supuesto de homocedasticidad y observando lo que nos arrojan las graficas podemos concluir que no se cumple este supuesto; ya que en los graficos a y b los datos no se destribuyen de manera equivalente, los datos tienden a estar en el lado izquierdo de las dos graficas.
car::qqPlot(modelo,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm')
## [1] 75 125
ks.test(modelo$residuals,pnorm,mean(modelo$residuals),sd(modelo$residuals))
## Warning in ks.test.default(modelo$residuals, pnorm, mean(modelo$residuals), :
## ties should not be present for the Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modelo$residuals
## D = 0.10619, p-value = 0.02101
## alternative hypothesis: two-sided
Con la grafica de normalidad, se puede observar que algunos de los datos no siguen el patron lineal; por lo cual se puede decir que no se cumple el supuesto de normalidad. Para confirmar lo observado en la grafica, se procede a realizar el test de Kolmogorov-Smirnov, en el cual se obtiene un valor p de 0.02101, el cual es menor al nivel de significancia y por ende se concluye que el modelo no cumple el supuesto de normalidad.
como el modelo no cumple ninguno de los supuestos, se procede a aplicar una tranformación al modelo:
boxcox = MASS::boxcox(modelo,lambda=seq(-3,3,length.out = 1000),
ylab='log-verosimilitud')
boxcox$x[boxcox$y ==max(boxcox$y)] # valor que maximiza la log-verosimilitud
## [1] 0.2012012
Una vez obtenido el \(\lambda = 0.2012012\), se procede a realizar la tranformación del modelo
modelonew<-lm(I(precio_millon^0.2012012)~Area_contruida,data = viviendanew)
summary(modelonew)
##
## Call:
## lm(formula = I(precio_millon^0.2012012) ~ Area_contruida, data = viviendanew)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57302 -0.10417 -0.00491 0.09291 0.42248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.032e+00 1.759e-02 172.31 <2e-16 ***
## Area_contruida 1.556e-03 7.002e-05 22.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1591 on 200 degrees of freedom
## Multiple R-squared: 0.7118, Adjusted R-squared: 0.7103
## F-statistic: 493.9 on 1 and 200 DF, p-value: < 2.2e-16
Cabe observar que con la transformación de box-cox realizada, el nuevo modelo sigue teniendo tanto el intercepto como la variable area construida como significativas; sin embargo el \(R^{2}\) ajustado disminuyo a un 71.03%.
Dicho lo anterior, se procede a realizar la validación de los supuestos para el nuevo modelo:
bptest(modelonew)
##
## studentized Breusch-Pagan test
##
## data: modelonew
## BP = 35.776, df = 1, p-value = 2.214e-09
Realizando el test de Breusch-pagan nuevamente, obtuvimos un valor p menor a 0.05 de tal manera que rechazamos la hipotesis nula y de esta manera no se cumple el supuesto de homocedasticidad.
###RESIDUOS
res.stud.base = studres(modelonew)
modelonew.fit.base= modelonew$fitted.values
par(mfrow=c(1,2))
plot(modelonew.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modelonew.fit.base), col = 2)
plot(modelonew.fit.base,abs(res.stud.base),
ylab='valor absoluto de los residuos estudentizados',
xlab='valores ajustados',main='(b)', col = "BLUE")
lines(lowess(abs(res.stud.base)~modelonew.fit.base), col = 2)
Comprobando el supuesto de homocedasticidad por medio de las graficas, se observa nuevamente una tendencia para el costado izquierdo evitando ser equivalente tanto en la grafica a como en la b. Con estó y el test anterior se confirma el que el supuesto de homocedasticidad no se logra cumplir en el nuevo modelo.
car::qqPlot(modelonew,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm')
## [1] 1 20
ks.test(modelonew$residuals,pnorm,mean(modelonew$residuals),sd(modelonew$residuals))
## Warning in ks.test.default(modelonew$residuals, pnorm,
## mean(modelonew$residuals), : ties should not be present for the Kolmogorov-
## Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modelonew$residuals
## D = 0.042428, p-value = 0.8603
## alternative hypothesis: two-sided
Observando la grafica de normalidad del nuevo modelo, se dice que los datos si presentan un patron lineal y por ende se concluye que se cumple el supuesto de normalidad. Por otro lado se confirma el supuesto al realizar el test de Kolmogorov-Smirnov, ya que se obtiene un valor p superior al nivel de significancia.
Se procede a realizar nuevamente la selección de variables para asi agregar la variable tipo de vivienda y lograr realizar una regresión categórica:
viviendanew2 <- vivienda %>%
filter(Barrio == "el ingenio")%>%
select(precio_millon,Area_contruida,Tipo,Barrio)
modelo2<-viviendanew2
modelo2<- dummy_cols(modelo2,
select_columns = "Tipo")
modelo2
## # A tibble: 202 × 6
## precio_millon Area_contruida Tipo Barrio Tipo_Apartamento Tipo_Casa
## <dbl> <dbl> <chr> <chr> <int> <int>
## 1 1450 1200 Casa el ingen… 0 1
## 2 290 100 Apartamento el ingen… 1 0
## 3 360 99 Apartamento el ingen… 1 0
## 4 550 197 Apartamento el ingen… 1 0
## 5 410 136 Apartamento el ingen… 1 0
## 6 390 198 Apartamento el ingen… 1 0
## 7 300 147 Apartamento el ingen… 1 0
## 8 250 97 Apartamento el ingen… 1 0
## 9 370 104 Apartamento el ingen… 1 0
## 10 1200 616 Casa el ingen… 0 1
## # … with 192 more rows
Obtenida la nueva base de datos con la variable categorica tipo de
vivienda y filtrada con el barrio “el ingenio”, Primero que todo se
procede a realizar una grafica descriptiva de la nueva base:
Con este grafico de cajas se puede concluir que las casas situadas en el barrio el ingenio, son mas costosas con respecto a los apartamentos situados en el mismo barrio.
Ahora, se procede a realizar el nuevo modelo, observando el comportamiento cuando el tipo de vivienda es “casa”:
modelocasa<-lm(precio_millon~Area_contruida+Tipo_Casa, data = modelo2)
summary(modelocasa)
##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Tipo_Casa, data = modelo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -277.45 -48.20 -9.33 53.92 491.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 196.11964 10.91567 17.967 < 2e-16 ***
## Area_contruida 1.02197 0.05616 18.196 < 2e-16 ***
## Tipo_Casa 99.24986 18.63125 5.327 2.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.68 on 199 degrees of freedom
## Multiple R-squared: 0.802, Adjusted R-squared: 0.8
## F-statistic: 403.1 on 2 and 199 DF, p-value: < 2.2e-16
Observando el modelo con el tipo de vivienda casa, se tiene que tanto el intercepto como la variable de area construida y el tipo casa son significativos para el modelo. Con respecto al coeficiente \(\beta_{1}\), se dice que por cada unidad de area construida va aumentar en 1.02197 para la variable respuesta (precio millon), y con respecto al coeficiente \(\beta_{2}\), se dice que si el tipo de vivienda es casa entonces tendra un incremento de 99.24986 en el precio de la vivienda. Ahora, observando el \(R^{2}\) del modelo, indica que el modelo calculado explica el 80.2% de la variabilidad presente en la variable respuesta, lo cual es algo muy bueno para el modelo y ademas, el \(R^{2}\) ajustado fue de 80%.
Se procede a realizar la validación de los supuestos:
##
## studentized Breusch-Pagan test
##
## data: modelocasa
## BP = 38.517, df = 2, p-value = 4.326e-09
Realizando el test de Breusch-pagan obtuvimos un valor p menor a 0.05 de tal manera que rechazamos la hipotesis nula y de esta manera no se cumple el supuesto de homocedasticidad.
###RESIDUOS
res.stud.base = studres(modelocasa)
modelocasa.fit.base= modelocasa$fitted.values
par(mfrow=c(1,2))
plot(modelocasa.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modelocasa.fit.base), col = 2)
plot(modelocasa.fit.base,abs(res.stud.base),
ylab='valor absoluto de los residuos estudentizados',
xlab='valores ajustados',main='(b)', col = "BLUE")
lines(lowess(abs(res.stud.base)~modelocasa.fit.base), col = 2)
Comprobando el supuesto de homocedasticidad y observando lo que nos arrojan las graficas podemos concluir que el modelo con tipo de vivienda casa en el barrio el ingenio no cumple este supuesto; ya que en los graficos a y b los datos no se distribuyen de manera equivalente, los datos tienden a estar en el lado izquierdo de las dos graficas.
car::qqPlot(modelocasa,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm')
## [1] 75 125
ks.test(modelocasa$residuals,pnorm,mean(modelocasa$residuals),sd(modelocasa$residuals))
## Warning in ks.test.default(modelocasa$residuals, pnorm,
## mean(modelocasa$residuals), : ties should not be present for the Kolmogorov-
## Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modelocasa$residuals
## D = 0.085829, p-value = 0.102
## alternative hypothesis: two-sided
Con la grafica de normalidad, se puede observar que la mayoria de los datos siguen el patron lineal; por lo cual se puede decir que se cumple el supuesto de normalidad. Para confirmar lo observado en la grafica, se procede a realizar el test de Kolmogorov-Smirnov, en el cual se obtiene un valor p de 0.102, el cual es mayor al nivel de significancia y por ende se concluye que el modelo con el tipo de vivienda casa cumple el supuesto de normalidad.
Ahora se procede a realizar el modelo para el tipo de vivienda apartamento y asi realizar la comparación:
modeloapart<-lm(precio_millon~Area_contruida+Tipo_Apartamento, data = modelo2)
summary(modeloapart)
##
## Call:
## lm(formula = precio_millon ~ Area_contruida + Tipo_Apartamento,
## data = modelo2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -277.45 -48.20 -9.33 53.92 491.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 295.36951 21.65885 13.637 < 2e-16 ***
## Area_contruida 1.02197 0.05616 18.196 < 2e-16 ***
## Tipo_Apartamento -99.24986 18.63125 -5.327 2.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.68 on 199 degrees of freedom
## Multiple R-squared: 0.802, Adjusted R-squared: 0.8
## F-statistic: 403.1 on 2 and 199 DF, p-value: < 2.2e-16
Con respecto al modelo con el tipo de vivienda apartamento, se tiene que tanto el intercepto como la variable de area construida y el tipo apartamento tambien son significativos para el modelo. Con respecto al coeficiente \(\beta_{1}\), se dice que por cada unidad de area construida va aumentar en 1.02197 para la variable respuesta(precio millon) al igual que con el modelo anterior, y con respecto al coeficiente \(\beta_{2}\), se dice que si el tipo de vivienda es aparmento entonces tendra un decrecimiento de 99.24986 en el precio de la vivienda (es la diferencia si el tipo de vivienda es casa o apartamento). Ahora, observando el \(R^{2}\) del modelo, indica que el modelo calculado al igual que el anterior, explica el 80.2% de la variabilidad presente en la variable respuesta, lo cual es algo muy bueno para el modelo y su \(R^{2}\) ajustado fue de 80%.
Se procede a realizar la validación de los supuestos:
##
## studentized Breusch-Pagan test
##
## data: modeloapart
## BP = 38.517, df = 2, p-value = 4.326e-09
Realizando el test de Breusch-pagan obtuvimos un valor p menor a 0.05 de tal manera que rechazamos la hipotesis nula y de esta manera no se cumple el supuesto de homocedasticidad.
###RESIDUOS
res.stud.base = studres(modeloapart)
modeloapart.fit.base= modeloapart$fitted.values
par(mfrow=c(1,2))
plot(modeloapart.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modeloapart.fit.base), col = 2)
plot(modeloapart.fit.base,abs(res.stud.base),
ylab='valor absoluto de los residuos estudentizados',
xlab='valores ajustados',main='(b)', col = "BLUE")
lines(lowess(abs(res.stud.base)~modeloapart.fit.base), col = 2)
Comprobando el supuesto de homocedasticidad y observando lo que nos arrojan las graficas podemos concluir que el modelo con tipo de vivienda apartamento en el barrio el ingenio no cumple este supuesto; ya que en los graficos a y b los datos no se distribuyen de manera equivalente, los datos tienden a estar en el lado izquierdo de las dos graficas al igual que el modelo de casa.
car::qqPlot(modeloapart,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm')
## [1] 75 125
ks.test(modeloapart$residuals,pnorm,mean(modeloapart$residuals),sd(modeloapart$residuals))
## Warning in ks.test.default(modeloapart$residuals, pnorm,
## mean(modeloapart$residuals), : ties should not be present for the Kolmogorov-
## Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modeloapart$residuals
## D = 0.085829, p-value = 0.102
## alternative hypothesis: two-sided
Con la grafica de normalidad, se puede observar que la mayoria de los datos siguen el patron lineal; por lo cual se puede decir que se cumple el supuesto de normalidad al igual que en modelo casa. Para confirmar lo observado en la grafica, se procede a realizar el test de Kolmogorov-Smirnov, en el cual se obtiene un valor p igual al modelo casa de 0.102, el cual es mayor al nivel de significancia y por ende se concluye que el modelo con el tipo de vivienda apartamento cumple el supuesto de normalidad.
Finalmente se concluye que los dos modelos, tanto el de casa como el de apartamento cumplen el supuesto de normalidad pero no cumplen el supuesto de homocedasticidad.
anova(modelocasa,modeloapart)
## Analysis of Variance Table
##
## Model 1: precio_millon ~ Area_contruida + Tipo_Casa
## Model 2: precio_millon ~ Area_contruida + Tipo_Apartamento
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 199 1937949
## 2 199 1937949 0 -2.3283e-10
Como se puede observar se obtiene un valor p inferior al nivel de significancia, por lo cual se puede decir que si se encuentra dierencias entre el precio en el tipo de vivienda casa y apartamento.
Comenzamos leyendo la base datosMe y vemos las primeras filas de la base:
## Masa Edad
## 1 106 43
## 2 106 41
## 3 97 47
## 4 113 46
## 5 96 45
## 6 119 41
Primero se ajusta el modelo cuadratico con la variable sin centrar
mod=lm(Masa ~ Edad+I(Edad^2), data = base)
summary(mod)
##
## Call:
## lm(formula = Masa ~ Edad + I(Edad^2), data = base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.086 -6.154 -1.088 6.220 20.578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 207.349608 29.225118 7.095 2.21e-09 ***
## Edad -2.964323 1.003031 -2.955 0.00453 **
## I(Edad^2) 0.014840 0.008357 1.776 0.08109 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.026 on 57 degrees of freedom
## Multiple R-squared: 0.7632, Adjusted R-squared: 0.7549
## F-statistic: 91.84 on 2 and 57 DF, p-value: < 2.2e-16
Vemos en los resultados del modelo que el intercepto y la variable Edad tienen un aporte significativo pero el termino de la Edad al cuadrado, no. Además de que el \(R^{2}\) ajustado fue de 75.49%.
Se realiza la comprobacion de los supuestos:
##
## studentized Breusch-Pagan test
##
## data: mod
## BP = 7.379, df = 2, p-value = 0.02498
Con el test de Breusch-Pagan se obtiene un valor p menor a 0.05, por lo que hay problemas de heterocedasticidad.
###RESIDUOS
res.stud.base = studres(mod)
mod.fit.base= mod$fitted.values
par(mfrow=c(1,2))
plot(mod.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~mod.fit.base), col = 2)
car::qqPlot(mod,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm',main='(b)')
## [1] 53 57
ks.test(mod$residuals,pnorm,mean(mod$residuals),sd(mod$residuals))
## Warning in ks.test.default(mod$residuals, pnorm, mean(mod$residuals),
## sd(mod$residuals)): ties should not be present for the Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: mod$residuals
## D = 0.085408, p-value = 0.7739
## alternative hypothesis: two-sided
En la figura anterior, en la grafica (a) de los residuos vemos que estos se distribuyen alrededor de cero y parece que unos valores ajustados pequeños, tienen residuos un poco altos. En la grafica (b) vemos el qqplot donde los puntos estan cerca de la linea, por lo que puede haber normalidad en los residuos. Tambien se realiza el test de Kolmogorov-Smirnov, en el cual se obtiene un valor p de 0.7739 y se concluye que se cumple el supuesto de normalidad.
Como el modelo no cumple el supuesto de normalidad, se realiza una transformacion por el metodo Box Cox:
boxcox=MASS::boxcox(mod,lambda=seq(-3,3,length.out = 1000),ylab='log-verosimilitud')
boxcox$x[boxcox$y ==max(boxcox$y)] # valor que maximiza la log-verosimilitud
## [1] 1.750751
Entonces se ajusta el modelo con la transformacion:
modT=lm(I(Masa^1.750751) ~ Edad+I(Edad^2), data = base)
summary(modT)
##
## Call:
## lm(formula = I(Masa^1.750751) ~ Edad + I(Edad^2), data = base)
##
## Residuals:
## Min 1Q Median 3Q Max
## -658.24 -304.50 -57.89 305.06 902.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9562.8811 1393.0831 6.865 5.35e-09 ***
## Edad -183.6967 47.8118 -3.842 0.000309 ***
## I(Edad^2) 1.0436 0.3983 2.620 0.011250 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 382.6 on 57 degrees of freedom
## Multiple R-squared: 0.7796, Adjusted R-squared: 0.7719
## F-statistic: 100.8 on 2 and 57 DF, p-value: < 2.2e-16
Vemos que el \(R^{2}\) ajustado fue de 77.19% y que todas las variables tienen un aporte significativo.
Se realiza la comprobacion de los supuestos:
##
## studentized Breusch-Pagan test
##
## data: modT
## BP = 2.425, df = 2, p-value = 0.2974
Con el test de Breusch-Pagan se obtiene un valor p mayor a 0.05, por lo que no hay problemas de heterocedasticidad.
###RESIDUOS
res.stud.base = studres(modT)
modT.fit.base= modT$fitted.values
par(mfrow=c(1,2))
plot(modT.fit.base,res.stud.base, ylab='residuos estudentizados',
xlab='valores ajustados',main='(a)', col = "BLUE")
abline(h=0,lty=2)
lines(lowess(res.stud.base~modT.fit.base), col = 2)
car::qqPlot(modT,xlab='cuantiles teóricos',ylab='residuos estudentizados', distribution = 'norm',main='(b)')
## [1] 44 53
ks.test(modT$residuals,pnorm,mean(modT$residuals),sd(modT$residuals))
## Warning in ks.test.default(modT$residuals, pnorm, mean(modT$residuals), : ties
## should not be present for the Kolmogorov-Smirnov test
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: modT$residuals
## D = 0.10061, p-value = 0.5781
## alternative hypothesis: two-sided
En la figura anterior, en la grafica (a) de los residuos vemos que estos se distribuyen alrededor de cero y no se nota un patron en los residuos. En la grafica (b) vemos el qqplot donde los puntos estan cerca de la linea, por lo que puede haber normalidad en los residuos. Tambien se realiza el test de Kolmogorov-Smirnov, en el cual se obtiene un valor p de 0.5781 y se concluye que se cumple el supuesto de normalidad.
vif(modT)
## Edad I(Edad^2)
## 128.256 128.256
Verificamos los VIF (Factor de Inflación de Varianza) y vemos que son muy altos, por lo que hay problemas de multicolinealidad.
Teniendo el modelo listo, ahora se centra la variable Edad restandole la media:
Edad.<-base$Edad - mean(base$Edad)
Se ajusta el modelo con la variable centrada
modC=lm(base$Masa ~ Edad.+I(Edad.^2))
summary(modC)
##
## Call:
## lm(formula = base$Masa ~ Edad. + I(Edad.^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.086 -6.154 -1.088 6.220 20.578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82.935749 1.543146 53.745 <2e-16 ***
## Edad. -1.183958 0.088633 -13.358 <2e-16 ***
## I(Edad.^2) 0.014840 0.008357 1.776 0.0811 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.026 on 57 degrees of freedom
## Multiple R-squared: 0.7632, Adjusted R-squared: 0.7549
## F-statistic: 91.84 on 2 and 57 DF, p-value: < 2.2e-16
Vemos que se consigue el mismo \(R^{2}\) y \(R^{2}\) ajustado teniendo en cuenta el modelo sin variables centradas y sin transformacion. Al comprobar los supuestos se obtuvieron los mismos resultados y la misma transformacion, por lo que se ajusta el modelo con la variable centrada y la transformacion, y este modelo cumpliria los supuestos.
modCT=lm(I(base$Masa^1.750751) ~ Edad.+I(Edad.^2))
summary(modCT)
##
## Call:
## lm(formula = I(base$Masa^1.750751) ~ Edad. + I(Edad.^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -658.24 -304.50 -57.89 305.06 902.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2299.1525 73.5576 31.26 <2e-16 ***
## Edad. -58.4949 4.2249 -13.85 <2e-16 ***
## I(Edad.^2) 1.0436 0.3983 2.62 0.0112 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 382.6 on 57 degrees of freedom
## Multiple R-squared: 0.7796, Adjusted R-squared: 0.7719
## F-statistic: 100.8 on 2 and 57 DF, p-value: < 2.2e-16
Ahora se realiza una grafica de este modelo ajustado, donde vemos que
la linea se ajusta a los puntos.
Analizando los coeficientes del modelo, tenemos que todos fueron significativos, ademas si tambien vemos la grafica anterior, la edad tiene un efecto negativo sobre la masa, pues a mayor edad, menor masa tendra la persona.
Se define un modelo reducido, el cual no cuenta con el termino cuadratico
modCTR=lm(I(base$Masa^1.750751) ~ Edad.)
Para probar si el termino cuadratico puede ser eliminado se plantean las siguientes hipotesis en terminos del coeficiente asociado al termino cuadratico: \(H_{o}: \beta_{2} = 0 \quad vs \quad H_{a}: \beta_{2} \neq 0\)
Esto lo probamos usando la funcion anova, donde ingresamos un modelo sin el término cuadrático, seguido del modelo que contiene el término cuadrático
anova(modCTR,modCT)
## Analysis of Variance Table
##
## Model 1: I(base$Masa^1.750751) ~ Edad.
## Model 2: I(base$Masa^1.750751) ~ Edad. + I(Edad.^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 9346446
## 2 57 8341859 1 1004587 6.8644 0.01125 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Como el valor p da menor a 0.05, se cuenta con evidencia para rechazar la hipotesis nula, por lo que el coeficiente \(\beta_{2}\) es diferente de cero. Por lo tanto, el término cuadrático no debería ser eliminado de este modelo.
Inicialmente vimos que los VIF del modelo sin la variable centrada eran mayores a 10 y esto denotaba problemas de multicolinealidad. Pero como se muestra a continuacion, al centrar la variable se consiguio eliminar el problema de multicolinealidad, pues los VIF son menores a 10.
vif(modCT)
## Edad. I(Edad.^2)
## 1.001473 1.001473
Tambien podemos ver las correlaciones, pues con la variable sin escalar se obtuvo una correlacion de 0.9960939.
cor(base$Edad,base$Edad^2)
## [1] 0.9960939
Pero cuando miramos la correlacion al centrar la variable, vemos que se reduce notoriamente, obteniendo -0.03835694.
cor(Edad.,Edad.^2)
## [1] -0.03835694
Finalmente, vimos que si se justifica el hecho de centrar la variable para eliminar la multicolinealidad entre dicha variable y su forma cuadrática.