Librerías

Subimos los siguientes paquetes:

library(plm)

## Warning: package 'plm' was built under R version 4.0.5

library(car)

## Loading required package: carData

library(ggstatsplot)

## Warning: package 'ggstatsplot' was built under R version 4.0.5

## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated

## Registered S3 methods overwritten by 'parameters':
##   method                           from      
##   as.double.parameters_kurtosis    datawizard
##   as.double.parameters_skewness    datawizard
##   as.double.parameters_smoothness  datawizard
##   as.numeric.parameters_kurtosis   datawizard
##   as.numeric.parameters_skewness   datawizard
##   as.numeric.parameters_smoothness datawizard
##   print.parameters_distribution    datawizard
##   print.parameters_kurtosis        datawizard
##   print.parameters_skewness        datawizard
##   summary.parameters_kurtosis      datawizard
##   summary.parameters_skewness      datawizard

## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

library(gplots)

## Warning: package 'gplots' was built under R version 4.0.5

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library(foreign)
library(lmtest)

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.0.5

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Ejemplo 1

Subimos a la sesión la base de datos: http://fmwww.bc.edu/ec-p/data/stockwatson/fatality.dta (noten que es un archivo .dta, es decir que fue creado en STATA, con lo que usamos la función read.dta del paquete foreign para leerlo)

imp_alc_accide_trans <- read.dta('http://fmwww.bc.edu/ec-p/data/stockwatson/fatality.dta')

Es una base de datos de accidentes de tráfico fatales en cada estado en EEUU excepto Alaska entre 1982 y 1988, que incluye datos sobre penas recibidas por las personas que fueron declaradas culpables de los mismos, la edad mínima legal para conducir, la tasa de desempleo, el pbi per cápita, e incluye la alícuota del impuesto sobre el consumo de cerveza. (Hay otras variables, pero nos concentraremos en estas).

El objetivo del trabajo fue determinar si el impuesto sobre el consumo de cerveza y la severidad de las penas influyen en la tasa de accidentes fatales. La base de datos (incluyendo las otra variables) fue estudiada en el paper: “Alcohol Policies and Highway Vehicle Fatalities” de C. J. Ruhm, Journal of Health Economics, 15(4), 1996, 435-54.

Las variables a utilizar son: - allmort = mortalidad por accidentes (en personas). Esta es nuestra variable de interés, pero la queremos expresar en TASA de mortalidad para poder compararla entre estados, con lo que la dividimos pop (población):

imp_alc_accide_trans$tasa_mort <- imp_alc_accide_trans$allmort/imp_alc_accide_trans$pop

state =estado
year =ejercicio
beertax=alícuota sobre cerveza (dólares por caja)
unrate=tasa de desempleo
perinc=pbi per cápita (usaremos el log)

imp_alc_accide_trans$logpbi <- log(imp_alc_accide_trans$perinc)

jaild=encarcelado
comserd=sentenciado a trabajo comunitario
vmiles=millas promedio recorridas por automovilistas
mlda=edad mínima legal para conducir (vamos a simplificarla a 18, 19, 20)

imp_alc_accide_trans$edadminima <- ifelse(imp_alc_accide_trans$mlda<19, 18,
                                          ifelse(imp_alc_accide_trans$mlda>=19 & imp_alc_accide_trans$mlda<20,
                                                 19,
                                                 ifelse(imp_alc_accide_trans$mlda>=20 & imp_alc_accide_trans$mlda<21,
                                                        20,21)))

Empecemos con una regresión simple entre accidentes mortales y la alícuota, en 2 ejercicios distintos, 1982 y 1988:

reg_imp_mort_82 <- lm(tasa_mort ~ beertax, 
                      data = imp_alc_accide_trans[imp_alc_accide_trans$year==1982,])
reg_imp_mort_88 <- lm(tasa_mort ~ beertax, 
                      data = imp_alc_accide_trans[imp_alc_accide_trans$year==1988,])

Podemos ver que en 1982 el impuesto no es significativo pero que con los datos de 1988 si lo es, pero POSITIVO!

Desde ya que una regresión PROMEDIA los coeficientes correspondientes a cada unidad (estado en este caso). Si es válido el supuesto habitual de que una variable influye sobre otra de similar manera para cada unidad de observación, entonces la regresión anterior es correcta. Pero si influye de manera diferente, entonces no lo es.

Si tenemos dos períodos, podemos calcular la diferencia en las tasas de mortalidad sobre la diferencia en las alícuotas PERO sabiendo que el termino de error va a ser la diferencia de los errores de las regresiónes para cada periodo. Esta regresión va a arrojar un estimador consistente.

imp_alc_accide_trans$dif_mort_88_82 <- imp_alc_accide_trans$tasa_mort[imp_alc_accide_trans$year==1988] - 
  imp_alc_accide_trans$tasa_mort[imp_alc_accide_trans$year==1982]
imp_alc_accide_trans$dif_imp_88_82  <- imp_alc_accide_trans$beertax[imp_alc_accide_trans$year==1988] - 
  imp_alc_accide_trans$beertax[imp_alc_accide_trans$year==1982]

reg_dif_88_82 <- lm(dif_mort_88_82 ~ dif_imp_88_82,
                    data = imp_alc_accide_trans)
summary(reg_dif_88_82)

## 
## Call:
## lm(formula = dif_mort_88_82 ~ dif_imp_88_82, data = imp_alc_accide_trans)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.227e-04 -9.619e-06  9.212e-06  2.229e-05  6.774e-05 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.204e-06  2.251e-06  -3.201   0.0015 ** 
## dif_imp_88_82 -1.041e-04  1.548e-05  -6.723 7.72e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.869e-05 on 334 degrees of freedom
## Multiple R-squared:  0.1192, Adjusted R-squared:  0.1166 
## F-statistic:  45.2 on 1 and 334 DF,  p-value: 7.719e-11

Vemos que un aumento en el impuesto REDUCE significativamente la tasa de mortalidad.

¿Y si T>2? Entonces tenemos datos en panel.

Primero veamos un modelo con efectos fijos:

reg_imp_mort_fe <-plm(tasa_mort ~ beertax, 
                      data=imp_alc_accide_trans,
                      model = 'within')
summary(reg_imp_mort_fe)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = tasa_mort ~ beertax, data = imp_alc_accide_trans, 
##     model = "within")
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Residuals:
##        Min.     1st Qu.      Median     3rd Qu.        Max. 
## -5.8696e-05 -8.2838e-06 -1.2701e-07  7.9545e-06  8.9780e-05 
## 
## Coefficients:
##            Estimate  Std. Error t-value Pr(>|t|)    
## beertax -6.5587e-05  1.8785e-05 -3.4915 0.000556 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.0785e-07
## Residual Sum of Squares: 1.0345e-07
## R-Squared:      0.040745
## Adj. R-Squared: -0.11969
## F-statistic: 12.1904 on 1 and 287 DF, p-value: 0.00055597

Vemos que a mayor impuesto menor tasa de accidentes

Hagamos el test de agrupamiento, para lo cual necesitamos correr el modelo de regresión sobre todo el periodo bajo estudio:

reg_imp_mort <- lm(tasa_mort ~ beertax, 
                   data = imp_alc_accide_trans)
pFtest(reg_imp_mort_fe, reg_imp_mort)

## 
##  F test for individual effects
## 
## data:  tasa_mort ~ beertax
## F = 52.179, df1 = 47, df2 = 287, p-value < 2.2e-16
## alternative hypothesis: significant effects

El p-value del test de agrupamiento es menor al 5%, con lo que rechazamos la hipótesis de agrupamiento y concluimos que la información del panel (tanto porque hay diferencias entre unidades o a través del tiempo o ambas) es importante. Debemos correr modelos de datos en panel.

Ya tenemos el modelo con efectos fijos. Corramos ahora un modelo con efectos aleatorios:

reg_imp_mort_re <- plm(tasa_mort ~ beertax, 
                       data=imp_alc_accide_trans, 
                       index=c("state", "year"), 
                       model="random")
summary(reg_imp_mort_re)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = tasa_mort ~ beertax, data = imp_alc_accide_trans, 
##     model = "random", index = c("state", "year"))
## 
## Balanced Panel: n = 48, T = 7, N = 336
## 
## Effects:
##                     var   std.dev share
## idiosyncratic 3.605e-10 1.899e-05 0.119
## individual    2.660e-09 5.158e-05 0.881
## theta: 0.8622
## 
## Residuals:
##        Min.     1st Qu.      Median     3rd Qu.        Max. 
## -4.7109e-05 -1.2005e-05 -2.1530e-06  9.1011e-06  9.6435e-05 
## 
## Coefficients:
##                Estimate  Std. Error z-value Pr(>|z|)    
## (Intercept)  2.0671e-04  9.9971e-06 20.6773   <2e-16 ***
## beertax     -5.2016e-06  1.2418e-05 -0.4189   0.6753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1.2648e-07
## Residual Sum of Squares: 1.2642e-07
## R-Squared:      0.00052508
## Adj. R-Squared: -0.0024674
## Chisq: 0.175467 on 1 DF, p-value: 0.6753

Con efectos aleatorios no encontramos significatividad en el coeficiente del impuesto.

Comparemos ambos modelos con el test de Hausman para decidirnos por uno u otro:

Test de Hausman

phtest(reg_imp_mort_fe, reg_imp_mort_re)

## 
##  Hausman Test
## 
## data:  tasa_mort ~ beertax
## chisq = 18.353, df = 1, p-value = 1.835e-05
## alternative hypothesis: one model is inconsistent

Con un p-value<5%, aceptamos el modelo de efectos FIJOS.

Podemos extraer las ordenadas al origen correspondientes a cada unidad de observación, es decir los efectos fijos por estado. Para ello se corre:

summary(fixef(reg_imp_mort_fe))

##      Estimate Std. Error t-value  Pr(>|t|)    
## AL 3.4776e-04 3.1336e-05 11.0980 < 2.2e-16 ***
## AZ 2.9099e-04 9.2539e-06 31.4452 < 2.2e-16 ***
## AR 2.8227e-04 1.3212e-05 21.3636 < 2.2e-16 ***
## CA 1.9682e-04 7.4007e-06 26.5943 < 2.2e-16 ***
## CO 1.9934e-04 8.0371e-06 24.8019 < 2.2e-16 ***
## CT 1.6154e-04 8.3913e-06 19.2506 < 2.2e-16 ***
## DE 2.1700e-04 7.7457e-06 28.0159 < 2.2e-16 ***
## FL 3.2095e-04 2.2151e-05 14.4890 < 2.2e-16 ***
## GA 4.0022e-04 4.6403e-05  8.6249 4.435e-16 ***
## ID 2.8086e-04 9.8767e-06 28.4368 < 2.2e-16 ***
## IL 1.5160e-04 7.8478e-06 19.3176 < 2.2e-16 ***
## IN 2.0161e-04 8.8672e-06 22.7364 < 2.2e-16 ***
## IA 1.9337e-04 1.0222e-05 18.9176 < 2.2e-16 ***
## KS 2.2544e-04 1.0863e-05 20.7528 < 2.2e-16 ***
## KY 2.2601e-04 8.0462e-06 28.0893 < 2.2e-16 ***
## LA 2.6305e-04 1.6266e-05 16.1714 < 2.2e-16 ***
## ME 2.3697e-04 1.6006e-05 14.8045 < 2.2e-16 ***
## MD 1.7712e-04 8.2458e-06 21.4800 < 2.2e-16 ***
## MA 1.3679e-04 8.6477e-06 15.8178 < 2.2e-16 ***
## MI 1.9931e-04 1.1663e-05 17.0888 < 2.2e-16 ***
## MN 1.5804e-04 9.3628e-06 16.8797 < 2.2e-16 ***
## MS 3.4485e-04 2.0936e-05 16.4717 < 2.2e-16 ***
## MO 2.1814e-04 9.2523e-06 23.5764 < 2.2e-16 ***
## MT 3.1172e-04 9.4413e-06 33.0170 < 2.2e-16 ***
## NE 1.9555e-04 1.0551e-05 18.5342 < 2.2e-16 ***
## NV 2.8769e-04 8.1056e-06 35.4922 < 2.2e-16 ***
## NH 2.2232e-04 1.4114e-05 15.7512 < 2.2e-16 ***
## NJ 1.3719e-04 7.3328e-06 18.7089 < 2.2e-16 ***
## NM 3.9040e-04 1.0154e-05 38.4492 < 2.2e-16 ***
## NY 1.2910e-04 7.5629e-06 17.0696 < 2.2e-16 ***
## NC 3.1872e-04 2.5173e-05 12.6609 < 2.2e-16 ***
## ND 1.8542e-04 1.0193e-05 18.1912 < 2.2e-16 ***
## OH 1.8032e-04 1.0193e-05 17.6910 < 2.2e-16 ***
## OK 2.9326e-04 1.8429e-05 15.9133 < 2.2e-16 ***
## OR 2.3096e-04 8.1175e-06 28.4526 < 2.2e-16 ***
## PA 1.7102e-04 8.6477e-06 19.7759 < 2.2e-16 ***
## RI 1.2126e-04 7.7533e-06 15.6395 < 2.2e-16 ***
## SC 4.0348e-04 3.5479e-05 11.3724 < 2.2e-16 ***
## SD 2.4739e-04 1.4121e-05 17.5195 < 2.2e-16 ***
## TN 2.6020e-04 9.1624e-06 28.3983 < 2.2e-16 ***
## TX 2.5602e-04 1.0853e-05 23.5889 < 2.2e-16 ***
## UT 2.3137e-04 1.5453e-05 14.9721 < 2.2e-16 ***
## VT 2.5116e-04 1.3973e-05 17.9751 < 2.2e-16 ***
## VA 2.1874e-04 1.4664e-05 14.9170 < 2.2e-16 ***
## WA 1.8181e-04 8.2328e-06 22.0836 < 2.2e-16 ***
## WV 2.5809e-04 1.0767e-05 23.9707 < 2.2e-16 ***
## WI 1.7184e-04 7.7457e-06 22.1848 < 2.2e-16 ***
## WY 3.2491e-04 7.2328e-06 44.9219 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

También es útil expresar los efectos fijos de cada unidad de observación en términos de los desvíos respectivos de la media de los efectos fijos. Para ello corremos:

summary(fixef(reg_imp_mort_fe, type = "dmean"))

##       Estimate  Std. Error  t-value  Pr(>|t|)    
## AL  1.1006e-04  3.1336e-05   3.5121 0.0005161 ***
## AZ  5.3283e-05  9.2539e-06   5.7579 2.188e-08 ***
## AR  4.4560e-05  1.3213e-05   3.3726 0.0008470 ***
## CA -4.0891e-05  7.4007e-06  -5.5254 7.380e-08 ***
## CO -3.8373e-05  8.0371e-06  -4.7744 2.876e-06 ***
## CT -7.6170e-05  8.3913e-06  -9.0773 < 2.2e-16 ***
## DE -2.0705e-05  7.7457e-06  -2.6731 0.0079465 ** 
## FL  8.3242e-05  2.2151e-05   3.7579 0.0002076 ***
## GA  1.6252e-04  4.6403e-05   3.5023 0.0005348 ***
## ID  4.3153e-05  9.8767e-06   4.3692 1.744e-05 ***
## IL -8.6107e-05  7.8478e-06 -10.9721 < 2.2e-16 ***
## IN -3.6099e-05  8.8672e-06  -4.0710 6.058e-05 ***
## IA -4.4338e-05  1.0222e-05  -4.3376 1.997e-05 ***
## KS -1.2266e-05  1.0863e-05  -1.1291 0.2597802    
## KY -1.1696e-05  8.0462e-06  -1.4536 0.1471402    
## LA  2.5344e-05  1.6266e-05   1.5581 0.1203234    
## ME -7.3922e-07  1.6006e-05  -0.0462 0.9631970    
## MD -6.0588e-05  8.2458e-06  -7.3478 2.102e-12 ***
## MA -1.0092e-04  8.6477e-06 -11.6700 < 2.2e-16 ***
## MI -3.8397e-05  1.1663e-05  -3.2922 0.0011185 ** 
## MN -7.9666e-05  9.3628e-06  -8.5087 9.901e-16 ***
## MS  1.0715e-04  2.0936e-05   5.1178 5.669e-07 ***
## MO -1.9571e-05  9.2523e-06  -2.1152 0.0352734 *  
## MT  7.4016e-05  9.4413e-06   7.8396 8.909e-14 ***
## NE -4.2162e-05  1.0551e-05  -3.9962 8.188e-05 ***
## NV  4.9978e-05  8.1056e-06   6.1659 2.373e-09 ***
## NH -1.5390e-05  1.4114e-05  -1.0904 0.2764610    
## NJ -1.0052e-04  7.3328e-06 -13.7083 < 2.2e-16 ***
## NM  1.5269e-04  1.0154e-05  15.0382 < 2.2e-16 ***
## NY -1.0861e-04  7.5629e-06 -14.3610 < 2.2e-16 ***
## NC  8.1009e-05  2.5173e-05   3.2180 0.0014386 ** 
## ND -5.2288e-05  1.0193e-05  -5.1299 5.345e-07 ***
## OH -5.7386e-05  1.0193e-05  -5.6301 4.288e-08 ***
## OK  5.5549e-05  1.8428e-05   3.0143 0.0028056 ** 
## OR -6.7445e-06  8.1175e-06  -0.8309 0.4067439    
## PA -6.6691e-05  8.6477e-06  -7.7120 2.049e-13 ***
## RI -1.1645e-04  7.7533e-06 -15.0194 < 2.2e-16 ***
## SC  1.6577e-04  3.5479e-05   4.6724 4.582e-06 ***
## SD  9.6834e-06  1.4121e-05   0.6858 0.4934227    
## TN  2.2490e-05  9.1624e-06   2.4546 0.0146996 *  
## TX  1.8308e-05  1.0853e-05   1.6869 0.0927111 .  
## UT -6.3395e-06  1.5453e-05  -0.4102 0.6819363    
## VT  1.3451e-05  1.3973e-05   0.9627 0.3365174    
## VA -1.8963e-05  1.4664e-05  -1.2931 0.1970014    
## WA -5.5897e-05  8.2328e-06  -6.7895 6.460e-11 ***
## WV  2.0380e-05  1.0767e-05   1.8929 0.0593813 .  
## WI -6.5871e-05  7.7457e-06  -8.5042 1.021e-15 ***
## WY  8.7205e-05  7.2328e-06  12.0568 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test de Breusch-Godfrey-Wooldridge correlación serial

pbgtest(reg_imp_mort_fe)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  tasa_mort ~ beertax
## chisq = 50.668, df = 7, p-value = 1.068e-08
## alternative hypothesis: serial correlation in idiosyncratic errors

Con p<0.05, HAY correlación parcial, pero no es preocupante en paneles cortos. De todas maneras conviene corregir las estimaciones. Para ello se corre un corrector llamado de Newey-West de la matriz de covarianzas. La correlación serial no altera los coeficientes obtenidos pero si los errores estándar y por ende el grado de significatividad de los coeficientes. En R corremos:

coeftest(reg_imp_mort_fe, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##            Estimate  Std. Error t value Pr(>|t|)  
## beertax -6.5587e-05  2.8837e-05 -2.2744  0.02368 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si comparamos este resultado con el del modelo con EF (beertax -6.5587e-05 1.8785e-05 -3.4915 0.000556 ***) vemos que el coeficiente es el mismo, pero cambia el error estándar y por ende el t-value y la prob. Aunque en este caso sigue siendo significativo. Si solamente hubiera correlación serial pero no dependencia transversal (es decir si solamente hubiera correlación entre períodos pero no entre unidades), este resultado corregido sería el que hay que informar y en el que habría que basarse para recomendaciones, conclusiones, etc. Pero tenemos que ver si es que no hay dependencia transversal.

Test de dependencia transversal

pcdtest(reg_imp_mort_fe)

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  tasa_mort ~ beertax
## z = 5.436, p-value = 5.449e-08
## alternative hypothesis: cross-sectional dependence

HAY dependencia transversal (p-value<0.05) Si solamente hubiera dependencia transversal pero no correlación serial (es decir si solamente hubiera correlación entre unidades pero no entre períodos), el siguiente sería el resultado a informar -vean que se agrupan por datos por ejercicio):

coeftest(reg_imp_mort_fe, vcov. = vcovHC, cluster = 'time')

## 
## t test of coefficients:
## 
##            Estimate  Std. Error t value  Pr(>|t|)    
## beertax -6.5587e-05  9.4573e-06 -6.9351 2.691e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pero en nuestro caso tenemos tanto DT como CS. Entonces hay que corregir la matriz de covarianzas por el método de Driscoll and Kraay, que controla la presencia de dependencia transversal y de correlación serial:

coeftest(reg_imp_mort_fe, vcov. = vcovSCC)

## 
## t test of coefficients:
## 
##            Estimate  Std. Error t value  Pr(>|t|)    
## beertax -6.5587e-05  1.0628e-05 -6.1712 2.303e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Este es el resultado a informar (además de los efectos fijos)

Ejemplo 2

Vamos a utilizar una base de datos que es una encuesta de opinión en 7 países entre 1990 y 1999 (la misma base que utilizaremos para ver Diferencias en Diferencias)

Panel <- read.dta("http://dss.princeton.edu/training/Panel101.dta")

Las variables son:

colnames(Panel)

## [1] "country" "year"    "y"       "y_bin"   "x1"      "x2"      "x3"     
## [8] "opinion" "op"

Simplemente corremos modelos de y sobre x1

1. Regresión lineal

ols <-lm(y ~ x1, data=Panel)
summary(ols)

## 
## Call:
## lm(formula = y ~ x1, data = Panel)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -9.546e+09 -1.578e+09  1.554e+08  1.422e+09  7.183e+09 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 1.524e+09  6.211e+08   2.454   0.0167 *
## x1          4.950e+08  7.789e+08   0.636   0.5272  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.028e+09 on 68 degrees of freedom
## Multiple R-squared:  0.005905,   Adjusted R-squared:  -0.008714 
## F-statistic: 0.4039 on 1 and 68 DF,  p-value: 0.5272

Vemos que x1 NO es significativa. Pero hay que ver si el supuesto de agrupamiento es aceptable

1. Efectos fijos

pan_fe <-plm(y ~ x1, data=Panel, model = 'within')
summary(pan_fe)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = y ~ x1, data = Panel, model = "within")
## 
## Balanced Panel: n = 7, T = 10, N = 70
## 
## Residuals:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -8.63e+09 -9.70e+08  5.40e+08  0.00e+00  1.39e+09  5.61e+09 
## 
## Coefficients:
##      Estimate Std. Error t-value Pr(>|t|)  
## x1 2475617827 1106675594   2.237  0.02889 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    5.2364e+20
## Residual Sum of Squares: 4.8454e+20
## R-Squared:      0.074684
## Adj. R-Squared: -0.029788
## F-statistic: 5.00411 on 1 and 62 DF, p-value: 0.028892

Cuando corremos el modelo con efectos fijos, resulta que x1 SI es significativa. Evidentemente, el modelo agregado no controlaba la heterogeneidad entre unidades de observación Podemos graficar las rectas de regresión de ambos modelos para comparar ambos modelos. Para ello, tenemos que correr el modelo con efectos fijos nuevamente, con lm() agregando una variable ficticia por unidad de observación, es decir, una ordenada por unidad. Fíjense que agregamos -1 para eliminar la ordenada al origen general, con lo que tendremos una ordenada por país.

fe_dum <-lm(y ~ x1 + factor(country) - 1, data=Panel)
Panel$yhat <- fe_dum$fitted  #valores ajustados de y

scatterplot(yhat~x1|country, data=Panel, boxplots=FALSE, 
            xlab="x1", ylab="yhat",smooth=FALSE)
abline(lm(y~x1, data = Panel),lwd=3, col="red")

Test de agrupamiento

pFtest(pan_fe, ols)

## 
##  F test for individual effects
## 
## data:  y ~ x1
## F = 2.9655, df1 = 6, df2 = 62, p-value = 0.01307
## alternative hypothesis: significant effects

Con p<0.5, se rechaza el agrupamiento. Es decir que se confirma que hay que correr un modelo de datos en panel. Corremos uno modelo con efectos aleatorios y lo comparamos con el de efectos fijos por medio del test de Hausman:

pan_re <- plm(y ~ x1, data=Panel, 
              index=c("country", "year"), model="random")
phtest(pan_fe, pan_re)

## 
##  Hausman Test
## 
## data:  y ~ x1
## chisq = 3.674, df = 1, p-value = 0.05527
## alternative hypothesis: one model is inconsistent

Con p>0.05, usar el modelo de efectos aleatorios.

Test de Breusch-Godfrey-Wooldridge correlación serial

pbgtest(pan_re)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  y ~ x1
## chisq = 10.274, df = 10, p-value = 0.4168
## alternative hypothesis: serial correlation in idiosyncratic errors

Con p>0.05, NO hay correlación serial.

Test de dependencia transversal

pcdtest(pan_re)

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  y ~ x1
## z = 1.5264, p-value = 0.1269
## alternative hypothesis: cross-sectional dependence

Con p>0.05, NO hay dependencia transversal.

Ejemplo 3

library(WDI)

WDI es un paquete para acceder a datos de los indicadores de desarrollo del Banco Mundial, el índice de desarrollo humano, indicadores de pobreza, etc. Variables: - NY.GDP.PCAP.KD - PBI per cápita (en dólares de 2015) - GB.XPD.RSDV.GD.ZS - Gasto en I&D (% PBI) - SH.XPD.GHED.GD.ZS - Gasto público en salud (% PBI) - SP.DYN.LE00.IN - Expectativa de vida al nacer

Queremos ver si la expectativa de vida al nacer depende del gasto en I&D, del gasto público en salud y del PBI per cápita entre el 2010 y el 2015 en Argentina, Brasil, Chile, México, Colombia, y Uruguay.

expectativa <- WDI(country=c('AR',
                             'BR',
                             'CL',
                             'CO',
                             'MX',
                             'UY'), start = 2010,
                   end=2015,
                   indicator=c('NY.GDP.PCAP.KD',
                               'GB.XPD.RSDV.GD.ZS',
                               'SH.XPD.GHED.GD.ZS',
                               'SP.DYN.LE00.IN'))
head(expectativa)

##   iso2c   country year NY.GDP.PCAP.KD GB.XPD.RSDV.GD.ZS SH.XPD.GHED.GD.ZS
## 1    AR Argentina 2010       13551.34           0.56104          5.568343
## 2    AR Argentina 2011       14200.27           0.56597          5.681173
## 3    AR Argentina 2012       13895.63           0.63491          6.134855
## 4    AR Argentina 2013       14071.51           0.61849          6.220423
## 5    AR Argentina 2014       13567.95           0.59396          6.327308
## 6    AR Argentina 2015       13789.06           0.62262          6.824373
##   SP.DYN.LE00.IN
## 1         75.721
## 2         76.124
## 3         76.467
## 4         76.491
## 5         76.755
## 6         76.760

colnames(expectativa)[4:7] <- c('pbi_pc',
                                'gasto_ID',
                                'gasto_salud',
                                'expectativa_nacer')
head(expectativa)

##   iso2c   country year   pbi_pc gasto_ID gasto_salud expectativa_nacer
## 1    AR Argentina 2010 13551.34  0.56104    5.568343            75.721
## 2    AR Argentina 2011 14200.27  0.56597    5.681173            76.124
## 3    AR Argentina 2012 13895.63  0.63491    6.134855            76.467
## 4    AR Argentina 2013 14071.51  0.61849    6.220423            76.491
## 5    AR Argentina 2014 13567.95  0.59396    6.327308            76.755
## 6    AR Argentina 2015 13789.06  0.62262    6.824373            76.760

Primero, corremos una regresión lineal:

expect_reglin <- lm(expectativa_nacer ~ log(pbi_pc) + 
               gasto_ID + 
               gasto_salud, 
             data=expectativa)
summary(expect_reglin)

## 
## Call:
## lm(formula = expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
##     data = expectativa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6614 -0.7565 -0.1456  0.7582  2.4300 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 52.02457    5.56551   9.348 1.15e-10 ***
## log(pbi_pc)  2.79155    0.61443   4.543 7.45e-05 ***
## gasto_ID    -3.05643    0.64431  -4.744 4.18e-05 ***
## gasto_salud -0.02715    0.18599  -0.146    0.885    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.222 on 32 degrees of freedom
## Multiple R-squared:  0.5928, Adjusted R-squared:  0.5546 
## F-statistic: 15.53 on 3 and 32 DF,  p-value: 2.074e-06

El pooling (que es como el anterior – POR QUE????):

expect_po <- plm(expectativa_nacer ~ log(pbi_pc) + 
                 gasto_ID + 
                 gasto_salud, 
               data=expectativa, 
               model = "pooling",
               index=c('country','year'))
summary(expect_po)

## Pooling Model
## 
## Call:
## plm(formula = expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
##     data = expectativa, model = "pooling", index = c("country", 
##         "year"))
## 
## Balanced Panel: n = 6, T = 6, N = 36
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -1.66144 -0.75650 -0.14556  0.75821  2.43000 
## 
## Coefficients:
##              Estimate Std. Error t-value  Pr(>|t|)    
## (Intercept) 52.024565   5.565509  9.3477 1.150e-10 ***
## log(pbi_pc)  2.791554   0.614427  4.5433 7.453e-05 ***
## gasto_ID    -3.056434   0.644313 -4.7437 4.179e-05 ***
## gasto_salud -0.027151   0.185985 -0.1460    0.8849    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    117.35
## Residual Sum of Squares: 47.789
## R-Squared:      0.59278
## Adj. R-Squared: 0.5546
## F-statistic: 15.5272 on 3 and 32 DF, p-value: 2.0736e-06

Test de Agrupamiento (‘poolability’)

pFtest(expectativa_nacer ~ log(pbi_pc) + 
           gasto_ID + 
           gasto_salud, 
         data=expectativa,
         model="within",
         index = c("country", 
                   "year"))

## 
##  F test for individual effects
## 
## data:  expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud
## F = 269.82, df1 = 5, df2 = 27, p-value < 2.2e-16
## alternative hypothesis: significant effects

Con p<0.5, se rechaza el agrupamiento.

Modelo 3 - FE (efectos fijos)

expectativa_fe <- plm(expectativa_nacer ~ log(pbi_pc) + 
                        gasto_ID + 
                        gasto_salud, 
                      data=expectativa,
                      model = "within",
                      index = c("country", 
                                "year"))
summary(expectativa_fe)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
##     data = expectativa, model = "within", index = c("country", 
##         "year"))
## 
## Balanced Panel: n = 6, T = 6, N = 36
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.3282930 -0.1321584 -0.0001745  0.1044142  0.3751160 
## 
## Coefficients:
##             Estimate Std. Error t-value  Pr(>|t|)    
## log(pbi_pc)  3.62129    0.84477  4.2867 0.0002065 ***
## gasto_ID     3.00224    0.73936  4.0606 0.0003769 ***
## gasto_salud  0.41289    0.14000  2.9492 0.0065054 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4.5217
## Residual Sum of Squares: 0.93765
## R-Squared:      0.79263
## Adj. R-Squared: 0.73119
## F-statistic: 34.4013 on 3 and 27 DF, p-value: 2.2845e-09

Modelo 4 - RE (efectos aleatorios; RE=random effects)

expectativa_re <- plm(expectativa_nacer ~ log(pbi_pc) + 
                        gasto_ID + 
                        gasto_salud, 
                      data=expectativa,
                      model = "random",
                      index = c("country", 
                                "year"))
summary(expectativa_re)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
##     data = expectativa, model = "random", index = c("country", 
##         "year"))
## 
## Balanced Panel: n = 6, T = 6, N = 36
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.03473 0.18635  0.01
## individual    3.46648 1.86185  0.99
## theta: 0.9592
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.467454 -0.097355 -0.012793  0.124402  0.358072 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 39.62756    7.47865  5.2988 1.166e-07 ***
## log(pbi_pc)  3.58213    0.83282  4.3012 1.699e-05 ***
## gasto_ID     2.52294    0.73517  3.4318 0.0005997 ***
## gasto_salud  0.44545    0.14232  3.1300 0.0017481 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4.7098
## Residual Sum of Squares: 1.2263
## R-Squared:      0.73963
## Adj. R-Squared: 0.71522
## Chisq: 90.9017 on 3 DF, p-value: < 2.22e-16

Test de Hausman

phtest(expectativa_fe, expectativa_re)

## 
##  Hausman Test
## 
## data:  expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud
## chisq = 43.934, df = 3, p-value = 1.558e-09
## alternative hypothesis: one model is inconsistent

Con p>0.05, se rechaza FE y se acepta RE

Test de Wooldridge correlación serial

pwartest(expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, data = expectativa, 
         index = c("country", "year"))

## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  plm.model
## F = 13.715, df1 = 1, df2 = 28, p-value = 0.0009254
## alternative hypothesis: serial correlation

Con p<0.05, se acepta que existe correlación serial.

Test de dependencia transversal

pcdtest(expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
        data = expectativa,
        index = c("country", 
                  "year"))

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud
## z = -0.64244, p-value = 0.5206
## alternative hypothesis: cross-sectional dependence

Con p>0.05, se acepta que no hay dependencia transversal.

Dado que hay correlación serial pero no dependencia transversal, corremos el corrector de la matriz de covarianza según el método de Newey-West:

expectativa_re_nw <- plm(expectativa_nacer ~ log(pbi_pc) + 
                        gasto_ID + 
                        gasto_salud, 
                      data=expectativa,
                      model = "random",
                      index = c("country", 
                                "year"),
                      vcov = vcovNW)
summary(expectativa_re_nw)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = expectativa_nacer ~ log(pbi_pc) + gasto_ID + gasto_salud, 
##     data = expectativa, model = "random", index = c("country", 
##         "year"), vcov = vcovNW)
## 
## Balanced Panel: n = 6, T = 6, N = 36
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.03473 0.18635  0.01
## individual    3.46648 1.86185  0.99
## theta: 0.9592
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.467454 -0.097355 -0.012793  0.124402  0.358072 
## 
## Coefficients:
##             Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 39.62756    7.47865  5.2988 1.166e-07 ***
## log(pbi_pc)  3.58213    0.83282  4.3012 1.699e-05 ***
## gasto_ID     2.52294    0.73517  3.4318 0.0005997 ***
## gasto_salud  0.44545    0.14232  3.1300 0.0017481 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    4.7098
## Residual Sum of Squares: 1.2263
## R-Squared:      0.73963
## Adj. R-Squared: 0.71522
## Chisq: 90.9017 on 3 DF, p-value: < 2.22e-16

Ejercicio

Corran el modelo del ejemplo 1 (el de los accidentes de tránsito) pero agregando las siguientes variables independientes a la tasa de impuesto:

unrate=tasa de desempleo
perinc=pbi per cápita (usaremos el log)
jaild=encarcelado
comserd=sentenciado a trabajo comunitario
vmiles=millas promedio recorridas por automovilistas
mlda=edad mínima legal para conducir

Analicen si se necesita considerar el carácter de panel de los datos, y si se lo necesita, corran un modelo con EF y uno con EA, compárenlos, elijan el que sugiera su test, determinen si hay correlación serial y dependencia transversal y apliquen el corrector de la matriz de covarianzas a partir de lo que sus resultados arrojen.

Expliquen los PASOS realizados y los RESULTADOS obtenidos en palabras.

Modelos para datos en panel

Jose Iparraguirre

2023-09-07

Librerías

Ejemplo 1

Ejemplo 2

Ejemplo 3

Ejercicio