#INTRODUCCIÓN
#Este trabajo se encuentra estructurado en tres partes, la primera parte trata de hacer un análisis descriptivo, después seleccionamos un modelo y lo estimamos, y por último, realizamos un adiscusión en comparación con el resultado que esperábamos, en conclusión, tenemos que intentar llegar a un modelo que sea lineal, homocedástico y normal.
#Teniendo en cuenta la variable explicativa, no tengo altas expectativas sobre que exista una relación causal, debido a que otras variables están más relacionadas con la variable explicada, por lo que probaré con las variables “Proporción de nacidos en el extranjero sobre la población total (Porcentaje)” y “Proporción de nacionales sobre la población total (Porcentaje)”, que parecen que está más relacionada con la variable explicada "Renta neta media anual por habitante (Euros)".
#ANÁLISIS DESCRIPTIVO
#Lo primero que tenemos que hacer es cargar las librerias que necesitamos.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(pxR)
## Loading required package: reshape2
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
## Loading required package: RJSONIO
## Loading required package: plyr
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:purrr':
##
## compact
library(printr)
## Registered S3 method overwritten by 'printr':
## method from
## knit_print.data.frame rmarkdown
library(mosaic)
## Loading required package: lattice
## Loading required package: ggformula
## Loading required package: ggstance
##
## Attaching package: 'ggstance'
## The following objects are masked from 'package:ggplot2':
##
## geom_errorbarh, GeomErrorbarh
##
## New to ggformula? Try the tutorials:
## learnr::run_tutorial("introduction", package = "ggformula")
## learnr::run_tutorial("refining", package = "ggformula")
## Loading required package: mosaicData
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Note: If you use the Matrix package, be sure to load it BEFORE loading mosaic.
##
## Attaching package: 'mosaic'
## The following object is masked from 'package:Matrix':
##
## mean
## The following object is masked from 'package:plyr':
##
## count
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
## The following object is masked from 'package:purrr':
##
## cross
## The following object is masked from 'package:ggplot2':
##
## stat
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median,
## prop.test, quantile, sd, t.test, var
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sum
library(mosaicModel)
## Loading required package: mosaicCore
##
## Attaching package: 'mosaicCore'
## The following object is masked from 'package:plyr':
##
## count
## The following objects are masked from 'package:dplyr':
##
## count, tally
## Loading required package: splines
##
## Attaching package: 'mosaicModel'
## The following objects are masked from 'package:mosaicCore':
##
## ci.mean, ci.median, ci.sd, coverage
library(broom)
library(ggfortify)
library(ggrepel)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:mosaicCore':
##
## logit
## The following objects are masked from 'package:mosaic':
##
## deltaMethod, logit
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
#Después, cargamos el archivo en formato zip que contiene el fichero pcAxis(.px) con los datos.
urban<-read.px("urbanaudit.px")%>%as.data.frame%>%as.tbl
#Una vez que tenemos los datos cargados, comenzamos a filtrar las variables, siendo la variable explicada "proporción de extranjeros sobre la población total (porcentaje) en el año 2015 para la población total" y la variable explicativa "renta neta media anual por habitante (euros) en el año 2015 para la población total"
urban1<-urban%>%filter(Sexo=="Total",Periodo=="2015",Indicadores%in%c("Proporción de extranjeros sobre la población total (Porcentaje)","Renta neta media anual por habitante (Euros)"))
#Para poder renombrar las variables, antes tenemos que pasar las variables de filas a columnas.
urban2<-urban1%>%spread(Indicadores,value)
#Ahora que las tenemos colocadas en columnas procedemos a renombrarlas con el comando "rename".
urban3<-urban2%>%dplyr::rename(Extr=`Proporción de extranjeros sobre la población total (Porcentaje)`)
urban4<-urban3%>%dplyr::rename(Renta=`Renta neta media anual por habitante (Euros)`)
#Hay que excluir el total nacional para no encontrarnos con datos ya repetidos.
urban5<-urban4[-c(1), ]
#Eliminamos los valores nulos del fichero.
urban6<-na.omit(urban5)
#Gráfico boxplot
#Podemos observar que nos encontramos con valores atípicos, resaltando un valor atípico más alto que el resto.
urban7<-urban6%>%ggplot(aes(x=Extr,y=Renta))+geom_boxplot()
urban7
## Warning: Continuous x aesthetic -- did you forget aes(group=...)?
#Función densidad
urban6%>%ggplot(aes(Renta))+geom_density()
#En este gráfico sobre la renta neta media anual por habitante en el año 2015, vemos que consta de una media unimodal y que es asimétrica hacia la derecha, lo que significa que la media aritmética es mayor que la mediana.
#Función de densidad
urban6%>%ggplot(aes(Extr))+geom_density()
#Ocurre algo muy parecido al gráfico de la variable explicativa, ya que también es unimodal, y asimétrica a la derecha, por lo que la media también es mayor a la mediana.
#Estadísticos descriptivos
summary(urban6)
|
2018: 0 |
Total :126 |
ES519C1 Albacete : 1 |
Min. : 1.150 |
Min. : 6753 |
|
2017: 0 |
Hombres: 0 |
41004 Alcalá de Guadaíra : 1 |
1st Qu.: 4.692 |
1st Qu.: 9510 |
|
2016: 0 |
Mujeres: 0 |
ES511C1 Alcalá de Henares: 1 |
Median : 8.745 |
Median :11090 |
|
2015:126 |
NA |
ES535C1 Alcobendas : 1 |
Mean :10.048 |
Mean :11262 |
|
NA |
NA |
ES517C1 Alcorcón : 1 |
3rd Qu.:12.672 |
3rd Qu.:12037 |
|
NA |
NA |
03009 Alcoy/Alcoi : 1 |
Max. :42.790 |
Max. :23861 |
|
NA |
NA |
(Other) :120 |
NA |
NA |
#Podemos comprobar como las conclusiones que sacamos de los gráficos boxplot son correctas, ya que nos indica que en ambas variables la media aritmética es superior a la mediana.
#SELECCIÓN DEL MODELO Y DIAGNÓSTICO
#Procedemos a realizar la estimación del modelo.
est1<-lm(Extr~Renta,urban6)
summary(est1)
##
## Call:
## lm(formula = Extr ~ Renta, data = urban6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.256 -5.247 -1.165 3.431 31.491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.4413851 2.7445593 4.897 2.96e-06 ***
## Renta -0.0003013 0.0002376 -1.268 0.207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.862 on 124 degrees of freedom
## Multiple R-squared: 0.01281, Adjusted R-squared: 0.004848
## F-statistic: 1.609 on 1 and 124 DF, p-value: 0.207
#Nos encontramos con que el R^2 ajustado del modelo es muy bajo (0.004848), lo cual quiere decir que la variable explicativa no es efectiva para explicar la variable explicada.
library(tigerstats)
## Loading required package: abd
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: grid
## Welcome to tigerstats!
## To learn more about this package, consult its website:
## http://homerhanumat.github.io/tigerstats
##
## Attaching package: 'tigerstats'
## The following object is masked from 'package:reshape2':
##
## tips
#Diagrama de dispersión
lmGC(Extr~Renta,data=urban6,graph=TRUE,check=TRUE)
##
## Linear Regression
##
## Correlation coefficient r = -0.1132
##
## Equation of Regression Line:
##
## Extr = 13.4414 + -3e-04 * Renta
##
## Residual Standard Error: s = 6.8616
## R^2 (unadjusted): R^2 = 0.0128
#Observamos que el modelo es lineal, ya que se encuentra dentro de las bandas de confianza.
#Gráficas para la valoración del modelo
autoplot(est1)
#En la gráfica "Residuals vs Fitted" vemos que la línea está muy cercana al cero, lo cual quiere decir que no hay evidencias de que no sea lineal.
#En segundo lugar, la gráfica "Scale-Location" nos hace pensar que el modelo es heterocedástico debido a que los puntos están muy juntos y conglomerados.
#Respecto a la gráfica Q-Q, nos informa sobre la normalidad del modelo, podemos decir que es normal ya que la línea se encuentra situada sobre la diagonal.
#Nos encontramos con algún atípico en la gráfica "Residuals vs Leverage".
#Mirando los gráficos, vemos que el modelo es lineal y heterocedástico, lo cual vamos a corroborarlo mediante el test Reset y el test de Breusch-Pagan.
#Linealidad
reset(est1)
##
## RESET test
##
## data: est1
## RESET = 1.709, df1 = 2, df2 = 122, p-value = 0.1853
#Como p-value= 0'1853>0'05 no rechazamos linalidad, por lo que el modelo es lineal.
#Homocedasticidad
bptest(est1)
##
## studentized Breusch-Pagan test
##
## data: est1
## BP = 11.233, df = 1, p-value = 0.0008036
#El modelo es heterocedástico debido a que p-value=0'0008036<0'05 rechazamos homocedasticidad, por lo que concluímos que el modelo es heterocedástico.
#Nos encontramos ante un modelo lineal y heterocedástico con un R^2 ajustado muy pequeño, por lo que voy a intentar mejorar el modelo añadiendo otra variable.
#Selección nuevo modelo
modelo1<-urban%>%filter(Sexo=="Total",Periodo==2015)
modelo2<-modelo1%>%filter(Indicadores=="Proporción de nacidos en el extranjero sobre la población total (Porcentaje)")
modelo3<-modelo2%>%spread(Indicadores,value)
modelo4<-modelo3%>%dplyr::rename(Nacext=`Proporción de nacidos en el extranjero sobre la población total (Porcentaje)`)
modelo5<-modelo4[-c(1), ]
modelo6<-na.omit(modelo5)
modelo7<-urban6%>%left_join(modelo6%>%select(Nacext,Nivel.territorial),by="Nivel.territorial")
#Gráfico boxplot
modelo8<-modelo7%>%ggplot(aes(x=Extr,y=Nacext))+geom_boxplot()
modelo8
## Warning: Continuous x aesthetic -- did you forget aes(group=...)?
#Función de densidad
modelo7%>%ggplot(aes(Nacext))+geom_density()
#Nos encontramos con un modelo unimodal y asimétrico a la derecha, por lo que la media sigue siendo mayor que la mediana.
#Estadísticos descriptivos
summary(modelo7)
|
2018: 0 |
Total :126 |
ES519C1 Albacete : 1 |
Min. : 1.150 |
Min. : 6753 |
Min. : 1.950 |
|
2017: 0 |
Hombres: 0 |
41004 Alcalá de Guadaíra : 1 |
1st Qu.: 4.692 |
1st Qu.: 9510 |
1st Qu.: 8.505 |
|
2016: 0 |
Mujeres: 0 |
ES511C1 Alcalá de Henares: 1 |
Median : 8.745 |
Median :11090 |
Median :11.865 |
|
2015:126 |
NA |
ES535C1 Alcobendas : 1 |
Mean :10.048 |
Mean :11262 |
Mean :13.494 |
|
NA |
NA |
ES517C1 Alcorcón : 1 |
3rd Qu.:12.672 |
3rd Qu.:12037 |
3rd Qu.:17.293 |
|
NA |
NA |
03009 Alcoy/Alcoi : 1 |
Max. :42.790 |
Max. :23861 |
Max. :47.010 |
|
NA |
NA |
(Other) :120 |
NA |
NA |
NA |
#Como hemos visto en la gráfica de la función de densidad de la nueva variable, afirmamos que la media es mayor que la mediana.
#Estimación del nuevo modelo
est3<-lm(Extr~Renta+Nacext,modelo7)
summary(est3)
##
## Call:
## lm(formula = Extr ~ Renta + Nacext, data = modelo7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0676 -0.7892 0.2715 0.8468 3.8152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.239e-01 7.209e-01 0.588 0.557599
## Renta -1.953e-04 5.719e-05 -3.415 0.000865 ***
## Nacext 8.762e-01 1.949e-02 44.945 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.651 on 123 degrees of freedom
## Multiple R-squared: 0.9433, Adjusted R-squared: 0.9424
## F-statistic: 1024 on 2 and 123 DF, p-value: < 2.2e-16
#Nos encontramos con un R^2 mucho mayor (0'9424), por lo que la variable explicativa es adecuada para explicar la variable adecuada.
library(ggfortify)
autoplot(est3)
#Si observamos la gráfica "Residuals vs Fitted", también nos afirma de, que no existen evidencias de que no sea lineal, ya que la línea continúa próxima al cero.
#En la gráfica "Scale-Location" vemos que el modelo es heterocedástico debido a que los puntos están muy juntos y conglomerados.
#En cuestión a la gráfica Q-Q, que nos informa sobre si el modelo es normal, vemos que sí lo es porque la línea QQ está ubicada encima de la diagonal.
#Nos encontramos con algún atípico en la gráfica "Residuals vs Leverage".
#Linealidad
reset(est3)
##
## RESET test
##
## data: est3
## RESET = 2.1141, df1 = 2, df2 = 121, p-value = 0.1252
#Como p-value= 0'1252>0'05 no rechazamos linealidad, por lo que el modelo es lineal.
#Homocedasticidad
bptest(est3)
##
## studentized Breusch-Pagan test
##
## data: est3
## BP = 16.986, df = 2, p-value = 0.0002049
#El modelo es heterocedástico debido a que p-value=0'0002049<0'05 rechazamos homocedasticidad.
#Este modelo es lineal y heterocedástico, pero con un R^2 ajustado mucho mayor, aún así, voy a intentar mejorar el modelo añadiendo otra variable para conseguir encontrar un modelo homocedástico.
#Selección nuevo modelo
tipo1<-urban%>%filter(Sexo=="Total",Periodo==2015)
tipo2<-tipo1%>%filter(Indicadores=="Proporción de nacionales sobre la población total (Porcentaje)")
tipo3<-tipo2%>%spread(Indicadores,value)
tipo4<-tipo3%>%dplyr::rename(Nac=`Proporción de nacionales sobre la población total (Porcentaje)`)
tipo5<-tipo4[-c(1), ]
tipo6<-na.omit(tipo5)
tipo7<-urban6%>%left_join(tipo6%>%select(Nac,Nivel.territorial),by="Nivel.territorial")
#Estimación modelo
est4<-lm(Extr~Renta+Nac,tipo7)
summary(est4)
##
## Call:
## lm(formula = Extr ~ Renta + Nac, data = tipo7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.945e-13 -2.557e-15 1.213e-15 5.257e-15 1.813e-14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.000e+02 2.189e-14 4.568e+15 <2e-16 ***
## Renta 4.305e-19 6.350e-19 6.780e-01 0.499
## Nac -1.000e+00 2.385e-16 -4.193e+15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.822e-14 on 123 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 8.904e+30 on 2 and 123 DF, p-value: < 2.2e-16
#El R^2 es muy bueno, por lo que la nueva variable explicativa elegida explica con efectividad la variable explicada.
library(ggfortify)
autoplot(est4)
#Si nos fijamos en la gráfica "Residuals vs Fitted" no encontramos evidencias de no linealidad.
#Respecto a la gráfica "Scale-Location" vemos que el modelo ya es homocedástico, debido a que no existen tendencias.
#Después, la gráfica Q-Q, nos información sobre que el modelo es normal del modelo, ya que la línea QQ está perfectamente situada sobre la diagonal.
#Si hacemos una visión global de los gráficos, nos encontramos con que el modelo es lineal y homocedástico. Vamos a comprobarlo con el test Reset y el test de Breush-Pagan.
#Linealidad
reset(est4)
##
## RESET test
##
## data: est4
## RESET = 0.21148, df1 = 2, df2 = 121, p-value = 0.8097
#Como p-value= 0'8097>0'05 no rechazamos la linealidad, por lo que continuamos concluyendo que el modelo es lineal.
#Homocedasticidad
bptest(est4)
##
## studentized Breusch-Pagan test
##
## data: est4
## BP = 0.66472, df = 2, p-value = 0.7172
#Vemos que este nuevo modelo es homocedástico ya que p-value=0'7172>0'05 no rechazamos homocedasticidad.
#DISCUSIÓN
#Como ya señalé en la introducción, buscamos un modelo que sea lineal y homocedástico, que la variable explicativa explicara efectivamente a la variable explicada, aunque, como era de esperar, la variable explicativa no lo hacía adecuadamente como hemos demostrado en el trabajo.
#Sin embargo, como intuí, las dos variables explicativas que he elegido explican de forma mucho más significativa la variable explicada, debido a que están más relacionadas.
#Como conclusión, he conseguido un modelo con variables más representativas, lineal, homocedástico y normal.