Actividad 1. Actividades R Parte 1 y 2

2024-02-15

EQUIPO 6

  • Alexa Mariana Marin Villar A00831342

  • Diego Martínez Ruibal A01740559

  • Oscar Emiliano Melendez Chavez A01276802

1.Obtener Indicadores del Banco Mundial

Librerías

#install.packages("WDI")
#install.packages("wbstats")
#install.packages ("gplots")
#install.packages("plm")
library(WDI)
library(wbstats)
library(tidyverse) ##limpia de datos incluye dplyr ytidyir
library(gplots)
library(plm)
library(mice)
#Imputacion PMM Mice
library(visdat)

Obtener información por país y eliminar columnas inservibles para el análisis

# PIB
# wb_data funcion para sacar la informacion del banco mundial los codigoss del pais son de dos letras
gdp_data <- wb_data(country = "MX", indicator = "NY.GDP.PCAP.CD", start_date = 2013, end_date = 2023)
gdp_data
## # A tibble: 10 × 9
##    iso2c iso3c country  date NY.GDP.PCAP.CD unit  obs_status footnote last_upd…¹
##    <chr> <chr> <chr>   <dbl>          <dbl> <chr> <chr>      <chr>    <date>    
##  1 MX    MEX   Mexico   2013         11317. <NA>  <NA>       <NA>     2023-12-19
##  2 MX    MEX   Mexico   2014         11490. <NA>  <NA>       <NA>     2023-12-19
##  3 MX    MEX   Mexico   2015         10098. <NA>  <NA>       <NA>     2023-12-19
##  4 MX    MEX   Mexico   2016          9153. <NA>  <NA>       <NA>     2023-12-19
##  5 MX    MEX   Mexico   2017          9693. <NA>  <NA>       <NA>     2023-12-19
##  6 MX    MEX   Mexico   2018         10130. <NA>  <NA>       <NA>     2023-12-19
##  7 MX    MEX   Mexico   2019         10435. <NA>  <NA>       <NA>     2023-12-19
##  8 MX    MEX   Mexico   2020          8895. <NA>  <NA>       <NA>     2023-12-19
##  9 MX    MEX   Mexico   2021         10359. <NA>  <NA>       <NA>     2023-12-19
## 10 MX    MEX   Mexico   2022         11497. <NA>  <NA>       <NA>     2023-12-19
## # … with abbreviated variable name ¹​last_updated
gdp_data<- select(gdp_data,-c(unit,obs_status,footnote))
gdp_data
## # A tibble: 10 × 6
##    iso2c iso3c country  date NY.GDP.PCAP.CD last_updated
##    <chr> <chr> <chr>   <dbl>          <dbl> <date>      
##  1 MX    MEX   Mexico   2013         11317. 2023-12-19  
##  2 MX    MEX   Mexico   2014         11490. 2023-12-19  
##  3 MX    MEX   Mexico   2015         10098. 2023-12-19  
##  4 MX    MEX   Mexico   2016          9153. 2023-12-19  
##  5 MX    MEX   Mexico   2017          9693. 2023-12-19  
##  6 MX    MEX   Mexico   2018         10130. 2023-12-19  
##  7 MX    MEX   Mexico   2019         10435. 2023-12-19  
##  8 MX    MEX   Mexico   2020          8895. 2023-12-19  
##  9 MX    MEX   Mexico   2021         10359. 2023-12-19  
## 10 MX    MEX   Mexico   2022         11497. 2023-12-19

Obtener información ditintos países

gdp_data2 <- wb_data(country = c("MX","EC","CA"), indicator = "NY.GDP.PCAP.CD", start_date = 2013, end_date = 2023)
gdp_data2
## # A tibble: 30 × 9
##    iso2c iso3c country  date NY.GDP.PCAP.CD unit  obs_status footnote last_upd…¹
##    <chr> <chr> <chr>   <dbl>          <dbl> <chr> <chr>      <chr>    <date>    
##  1 CA    CAN   Canada   2013         52635. <NA>  <NA>       <NA>     2023-12-19
##  2 CA    CAN   Canada   2014         50956. <NA>  <NA>       <NA>     2023-12-19
##  3 CA    CAN   Canada   2015         43596. <NA>  <NA>       <NA>     2023-12-19
##  4 CA    CAN   Canada   2016         42316. <NA>  <NA>       <NA>     2023-12-19
##  5 CA    CAN   Canada   2017         45129. <NA>  <NA>       <NA>     2023-12-19
##  6 CA    CAN   Canada   2018         46549. <NA>  <NA>       <NA>     2023-12-19
##  7 CA    CAN   Canada   2019         46374. <NA>  <NA>       <NA>     2023-12-19
##  8 CA    CAN   Canada   2020         43350. <NA>  <NA>       <NA>     2023-12-19
##  9 CA    CAN   Canada   2021         52359. <NA>  <NA>       <NA>     2023-12-19
## 10 CA    CAN   Canada   2022         54918. <NA>  <NA>       <NA>     2023-12-19
## # … with 20 more rows, and abbreviated variable name ¹​last_updated
gdp_date<- select(gdp_data2,-c(unit,obs_status,footnote))
gdp_date
## # A tibble: 30 × 6
##    iso2c iso3c country  date NY.GDP.PCAP.CD last_updated
##    <chr> <chr> <chr>   <dbl>          <dbl> <date>      
##  1 CA    CAN   Canada   2013         52635. 2023-12-19  
##  2 CA    CAN   Canada   2014         50956. 2023-12-19  
##  3 CA    CAN   Canada   2015         43596. 2023-12-19  
##  4 CA    CAN   Canada   2016         42316. 2023-12-19  
##  5 CA    CAN   Canada   2017         45129. 2023-12-19  
##  6 CA    CAN   Canada   2018         46549. 2023-12-19  
##  7 CA    CAN   Canada   2019         46374. 2023-12-19  
##  8 CA    CAN   Canada   2020         43350. 2023-12-19  
##  9 CA    CAN   Canada   2021         52359. 2023-12-19  
## 10 CA    CAN   Canada   2022         54918. 2023-12-19  
## # … with 20 more rows

Generar un conjunto de datos de panel solo con las variables necesarias

panel <- select(gdp_date,country,date,NY.GDP.PCAP.CD)
panel
## # A tibble: 30 × 3
##    country  date NY.GDP.PCAP.CD
##    <chr>   <dbl>          <dbl>
##  1 Canada   2013         52635.
##  2 Canada   2014         50956.
##  3 Canada   2015         43596.
##  4 Canada   2016         42316.
##  5 Canada   2017         45129.
##  6 Canada   2018         46549.
##  7 Canada   2019         46374.
##  8 Canada   2020         43350.
##  9 Canada   2021         52359.
## 10 Canada   2022         54918.
## # … with 20 more rows

1.1 Conjunto de Datos de Panel con Indicadores del Banco Mundial

Elabora un conjunto de datos de panel utilizando R. Elige una variable endógena dependiendo del tema que sea por equipo y selecciona otras variables como factores.

Género

Variable principal:

  1. GENDER: SL.IND.EMPL.FE.ZS = Employment in industry, female (% of female employment) (modeled ILO estimate)

Factores

  1. ECONOMY = NY.GDP.PCAP.CD = GDP per capita (current US)

  2. EDUCATIONSE.XPD.TOTL.GD.ZS = Government expenditure on education, total (% of GDP)

  3. EDUCATION SL.TLF.TOTL.FE.ZS = Labor force, female (% of total labor force)

  4. SOCIAL DEVELOPMENT SP.DYN.LE00.FE.IN = Life expectancy at birth, female (years)

  5. SOCIAL DEVELOPMENT SG.GEN.PARL.ZS= Proportion of seats held by women in national parliaments (%)

exercise_1 <- wb_data(country = c("MX","IT","JP","BR"), indicator = c("SL.IND.EMPL.FE.ZS","NY.GDP.PCAP.CD","SE.XPD.TOTL.GD.ZS","SL.TLF.TOTL.FE.ZS","SP.DYN.LE00.FE.IN","SG.GEN.PARL.ZS"), start_date = 2011, end_date = 2021)
exercise_1
## # A tibble: 44 × 10
##    iso2c iso3c country  date NY.GDP.PC…¹ SE.XP…² SG.GE…³ SL.IN…⁴ SL.TL…⁵ SP.DY…⁶
##    <chr> <chr> <chr>   <dbl>       <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 BR    BRA   Brazil   2011      13201.    5.74    8.58    12.1    42.2    76.7
##  2 BR    BRA   Brazil   2012      12328.    5.86    8.58    12.8    42.1    76.9
##  3 BR    BRA   Brazil   2013      12259.    5.84    8.58    12.4    42.1    77.1
##  4 BR    BRA   Brazil   2014      12071.    5.95    9.94    12.1    42.2    77.4
##  5 BR    BRA   Brazil   2015       8783.    6.24    9.94    11.8    42.6    77.6
##  6 BR    BRA   Brazil   2016       8681.    6.31    9.94    10.9    42.7    77.8
##  7 BR    BRA   Brazil   2017       9897.    6.32   10.7     10.8    43.2    78.0
##  8 BR    BRA   Brazil   2018       9121.    6.09   15.0     10.8    43.4    78.3
##  9 BR    BRA   Brazil   2019       8845.    5.96   15.0     10.7    43.7    78.5
## 10 BR    BRA   Brazil   2020       6924.    5.77   14.6     10.7    42.6    77.4
## # … with 34 more rows, and abbreviated variable names ¹​NY.GDP.PCAP.CD,
## #   ²​SE.XPD.TOTL.GD.ZS, ³​SG.GEN.PARL.ZS, ⁴​SL.IND.EMPL.FE.ZS,
## #   ⁵​SL.TLF.TOTL.FE.ZS, ⁶​SP.DYN.LE00.FE.IN

Imputar valores nulos

exercise_1 <- mice(exercise_1, m=2, maxit=1, meth='cart', seed=500)
## 
##  iter imp variable
##   1   1  SE.XPD.TOTL.GD.ZS
##   1   2  SE.XPD.TOTL.GD.ZS
exercise_1 <- complete(exercise_1)
vis_miss(exercise_1)

1.2 Gráficas de Heterogeneidad 

Observación de Heterogeneidad entre paises

# Variable de interes y por año
plotmeans(SL.IND.EMPL.FE.ZS ~ country, data= exercise_1, main= " Heterogeneidad entre países")

Observación de Heterogeneidad entre el paso de los años

plotmeans(SL.IND.EMPL.FE.ZS ~ date, data= exercise_1, main= " Heterogeneidad entre los años")

Heterogeneidad = variabilidad de los datos Homegeneidad = no hay varaqcion en los datos

1.3 Preguntas

1.¿Las lineas que une los promedios es horizontal o tiene muchos picos?

Las líneas presentan algunas variacione en el caso de los años, pero ninguna de ellas es muy extrema. Sin embargo, se puede observar que existe una mayor variación a partir del año 2015. En cuanto a los países, se presenta una diferencia notable en los promedios del indicador, ya que México es el país con la mayor tasa de empleo de mujeres. Este resultado resulta sorprendente y puede llevar a la formulación de diversas hipótesis sobre lo que podría estar ocurriendo en cada país.

2.¿Los intervalos de confianza miden o mismo, o estan desfasados?

En el caso de los países, existe un desfase con respecto a un país en particular. Por otro lado, en el caso de los años, no se observa una variación tan amplia, pero sí se aprecia cierta variabilidad a lo largo del tiempo.

3.¿Investiga elconcepto de heterogeneidad y determina si lo que se ve en las graficas es deseable o no deseable?

De acuerdo con la investigación, se encontró que se busca la heterogeneidad en los datos para permitir una exploración más profunda, comparación de datos y descubrimientos de ellos. En este caso, el porcentaje de empleo de mujeres en todo el mundo puede indicar las diferencias culturakles, sociales y economcias que existen entre países y a lo largo de los años, facilitando así una comprensión más amplia del comportamiento de la variable.

Investigar estas diferencias/ variabilidad de los datos es crucial para comprender por qué la variable y los factores pueden comportarse de manera diferente en diferentes situaciones. Puede revelar factores específico como características demográficas, condiciones socioeconómicas o diferencias culturales y del resultado que se puede esperar la próxima vez que se aplique.

2.Prueba de cambiar los años agrupados para observar una mayor variación

panel2 <- select(exercise_1,country,date,SL.IND.EMPL.FE.ZS,NY.GDP.PCAP.CD,SE.XPD.TOTL.GD.ZS,SL.TLF.TOTL.FE.ZS,SP.DYN.LE00.FE.IN,SG.GEN.PARL.ZS)
panel_prueba <- subset(panel2,date == 2011 | date == 2014 | date == 2018| date == 2021)
View(panel2)

2.1 Plot para observar los valores agrupados

plotmeans(SL.IND.EMPL.FE.ZS ~ date, data= panel_prueba, main= " Heterogeneidad entre los años")

3. Modelos con Indicadores del Banco Mundial

En este modelado de datos se decidio utilizar todos los años (2011-2021).

Modelo 1. Regresion agrupada (pooled)

modelo1 <- plm(SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN , data = panel2, model = "pooling")
summary(modelo1)
## Pooling Model
## 
## Call:
## plm(formula = SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + 
##     SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, data = panel2, model = "pooling")
## 
## Balanced Panel: n = 4, T = 11, N = 44
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.20946 -0.21661  0.02643  0.32160  0.72448 
## 
## Coefficients:
##                      Estimate  Std. Error  t-value  Pr(>|t|)    
## (Intercept)        6.7708e+01  4.7447e+00  14.2702 < 2.2e-16 ***
## NY.GDP.PCAP.CD     2.3459e-05  1.9073e-05   1.2300  0.226080    
## SE.XPD.TOTL.GD.ZS -1.6998e+00  1.8753e-01  -9.0642 3.831e-11 ***
## SL.TLF.TOTL.FE.ZS -6.4823e-01  5.2196e-02 -12.4192 3.953e-15 ***
## SP.DYN.LE00.FE.IN -2.3992e-01  6.9616e-02  -3.4463  0.001375 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    186.14
## Residual Sum of Squares: 8.2353
## R-Squared:      0.95576
## Adj. R-Squared: 0.95122
## F-statistic: 210.624 on 4 and 39 DF, p-value: < 2.22e-16

Modelo 2. Efectos fijos (within)

modelo1.2 <- plm(SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN , data = panel2, model = "within")
summary(modelo1.2)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + 
##     SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, data = panel2, model = "within")
## 
## Balanced Panel: n = 4, T = 11, N = 44
## 
## Residuals:
##        Min.     1st Qu.      Median     3rd Qu.        Max. 
## -0.89921750 -0.12122418  0.00099016  0.29311858  0.70856979 
## 
## Coefficients:
##                      Estimate  Std. Error t-value  Pr(>|t|)    
## NY.GDP.PCAP.CD     5.4739e-05  2.6986e-05  2.0285 0.0499613 *  
## SE.XPD.TOTL.GD.ZS -1.0574e+00  2.8566e-01 -3.7015 0.0007132 ***
## SL.TLF.TOTL.FE.ZS -3.4940e-01  1.2066e-01 -2.8957 0.0063942 ** 
## SP.DYN.LE00.FE.IN -2.1310e-01  8.4392e-02 -2.5251 0.0161108 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    12.893
## Residual Sum of Squares: 5.7541
## R-Squared:      0.55369
## Adj. R-Squared: 0.46691
## F-statistic: 11.1654 on 4 and 36 DF, p-value: 5.4139e-06

Prueba pF

#Efecto fijos primero y pooling despues
#significant effects significa que si hay diferencia entre los modelos (uno es mejor que otro) non significant puedes usar el que sea de los modelos

pFtest(modelo1.2,modelo1)
## 
##  F test for individual effects
## 
## data:  SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS +  ...
## F = 5.1746, df1 = 3, df2 = 36, p-value = 0.004469
## alternative hypothesis: significant effects

En este caso, existen diferencias entre ambos modelos fijos, por lo cual el modelo pooled es el elegido, ya que tiene una R cuadrada ajustada con un resultado de 0.95122.

Modelo 3. Efectos aleatorios (random) - Metodo walhus

modelo1.3 <- plm(SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, 
                 data = panel2,
                 model = "random", 
                 random.method = "walhus")
summary(modelo1.3)
## Oneway (individual) effect Random Effect Model 
##    (Wallace-Hussain's transformation)
## 
## Call:
## plm(formula = SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + 
##     SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, data = panel2, model = "random", 
##     random.method = "walhus")
## 
## Balanced Panel: n = 4, T = 11, N = 44
## 
## Effects:
##                    var  std.dev share
## idiosyncratic 0.186514 0.431873 0.997
## individual    0.000651 0.025515 0.003
## theta: 0.01866
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -1.209275 -0.219136  0.025578  0.320579  0.722151 
## 
## Coefficients:
##                      Estimate  Std. Error  z-value  Pr(>|z|)    
## (Intercept)        6.7590e+01  4.7491e+00  14.2320 < 2.2e-16 ***
## NY.GDP.PCAP.CD     2.3469e-05  1.9046e-05   1.2322 0.2178672    
## SE.XPD.TOTL.GD.ZS -1.6934e+00  1.8804e-01  -9.0055 < 2.2e-16 ***
## SL.TLF.TOTL.FE.ZS -6.4897e-01  5.2506e-02 -12.3597 < 2.2e-16 ***
## SP.DYN.LE00.FE.IN -2.3845e-01  6.9650e-02  -3.4236 0.0006181 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    179.73
## Residual Sum of Squares: 8.2063
## R-Squared:      0.95434
## Adj. R-Squared: 0.94966
## Chisq: 815.161 on 4 DF, p-value: < 2.22e-16

Modelo 3.2 Efectos aleatorios (random) - Metodo ameniya

#modelo1.4 <- plm(SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN,  data = panel2, model = "random", random.method = "ameniya")
#summary(modelo1.4)

#El modelo ameniya marca un error el cual es "Error in ercomp.formula(object, data, effect = effect, method = method,  :ameniya is not a relevant method"

El modelo Amemiya tiene un error, el cual es “Error in ercomp.formula(object, data, effect = effect, method = method, : ameniya is not a relevant method” por lo cual, en este caso, se decidió omitir el modelo.

Modelo 3.3

modelo1.5 <- plm(SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD 
                 + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, 
                 data = panel2,
                 model = "random", 
                 random.method = "nerlove")
summary(modelo1.5)
## Oneway (individual) effect Random Effect Model 
##    (Nerlove's transformation)
## 
## Call:
## plm(formula = SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + 
##     SL.TLF.TOTL.FE.ZS + SP.DYN.LE00.FE.IN, data = panel2, model = "random", 
##     random.method = "nerlove")
## 
## Balanced Panel: n = 4, T = 11, N = 44
## 
## Effects:
##                  var std.dev share
## idiosyncratic 0.1308  0.3616 0.097
## individual    1.2243  1.1065 0.903
## theta: 0.9019
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -1.028068 -0.175332  0.012991  0.257813  0.641837 
## 
## Coefficients:
##                      Estimate  Std. Error z-value  Pr(>|z|)    
## (Intercept)        5.3323e+01  6.5836e+00  8.0993 5.525e-16 ***
## NY.GDP.PCAP.CD     4.4400e-05  2.1832e-05  2.0337 0.0419825 *  
## SE.XPD.TOTL.GD.ZS -1.1016e+00  2.5386e-01 -4.3393 1.430e-05 ***
## SL.TLF.TOTL.FE.ZS -4.1467e-01  1.0796e-01 -3.8409 0.0001226 ***
## SP.DYN.LE00.FE.IN -2.2213e-01  7.3811e-02 -3.0095 0.0026168 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    14.559
## Residual Sum of Squares: 6.052
## R-Squared:      0.5843
## Adj. R-Squared: 0.54166
## Chisq: 54.8177 on 4 DF, p-value: 3.5477e-11

Prueba de Hausman

En este caso, el Modelo 1 del tipo fijo pooled cuenta con un coeficiente de determinación (R cuadrada) de 0.95122, mientras que el Modelo 1.3 del tipo aleatorio walhus tiene un valor de 0.94966. Aunque los resultados son muy parecidos, el test de Hausman nos permite observar que hay un mejor modelo y en este caso es el del tipo Pooled.

#El mejor de los tres aleatoruos que haya salido primero en comparativa con el mejor agrupado o fijo.
#hausman = fijos y aleatorios
phtest(modelo1.3,modelo1)
## 
##  Hausman Test
## 
## data:  SL.IND.EMPL.FE.ZS ~ NY.GDP.PCAP.CD + SE.XPD.TOTL.GD.ZS + SL.TLF.TOTL.FE.ZS +  ...
## chisq = 0.24411, df = 4, p-value = 0.9931
## alternative hypothesis: one model is inconsistent