Comparación de grupos en modelos FE de medidas repetidas, con un componente de error anidado, estimado utilizando plm

8

He estimado algunas medidas repetidas de modelos de efectos fijos, con un componente de error anidado, basado en variables de agrupación, es decir, modelos no anidados, utilizando . Ahora estoy interesado en

  1. pruebe si los modelos completos son significativamente diferentes, es decir, donde es el modelo completo para y es el modelo completo para y
    Ho:βFmimetrounalmi=βMETROunalmi
    βFmimetrounalmiFemalesβMETROunalmiMales
  2. posteriormente pruebe los coeficientes de regresión seleccionados entre dos grupos, es decir, donde es el coeficiente de regresión para las mujeres at , y es el coeficiente de regresión para hombres en .
    Ho:βFmimetrounalmi==ymiunar1,5=βMETROunalmi==ymiunar1,5
    βFmimetrounalmi==ymiunar1,5year1.5βMETROunalmi==ymiunar1,5year1.5

Ilustraré la situación usando el siguiente ejemplo de trabajo,

Primero, se necesitan algunos paquetes,

# install.packages(c("plm","texreg","tidyverse","lmtest"), dependencies = TRUE)
library(plm); library(lmtest); require(tidyverse)

Segundo, algo de preparación de datos,

data(egsingle, package = "mlmRev")
dta <-  egsingle %>% mutate(Female = recode(female,.default = 0L,`Female` = 1L))

Tercero, calculo un conjunto de modelos para cada género en los datos

MoSpc <- as.formula(math ~ Female + size + year)
dfMo = dta %>% group_by(female) %>%
    do(fitMo = plm(update(MoSpc, . ~ . -Female), 
       data = ., index = c("childid", "year", "schoolid"), model="within") )

Adelante, veamos los dos modelos estimados,

texreg::screenreg(dfMo[[2]], custom.model.names = paste0('FE: ', dfMo[[1]]))
#> ===================================
#>            FE: Female   FE: Male   
#> -----------------------------------
#> year-1.5      0.79 ***     0.88 ***
#>              (0.07)       (0.10)   
#> year-0.5      1.80 ***     1.88 ***
#>              (0.07)       (0.10)   
#> year0.5       2.51 ***     2.56 ***
#>              (0.08)       (0.10)   
#> year1.5       3.04 ***     3.17 ***
#>              (0.08)       (0.10)   
#> year2.5       3.84 ***     3.98 ***
#>              (0.08)       (0.10)   
#> -----------------------------------
#> R^2           0.77         0.79    
#> Adj. R^2      0.70         0.72    
#> Num. obs.  3545         3685       
#> ===================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05    #> 

Ahora, quiero probar si estos dos modelos (OLS lineales) son significativamente diferentes, cf. punto 1 arriba. Miré a mi alrededor SO e Internet y algunos sugieren que necesito usar plm::pFtest(), también sugerí aquí , que he intentado, pero no estoy convencido. Me hubiera imaginado alguna prueba para modelos no anidados, posible prueba de Cox lmtest::coxtest, pero no estoy seguro en absoluto. Si alguien aquí podría ayudarme.

Lo intenté,

plm::pFtest(dfMo[[1,2]], dfMo[[2,2]])
# >
# > F test for individual effects
# >
# >data:  update(MoSpc, . ~ . - Female)
# >F = -0.30494, df1 = 113, df2 = 2693, p-value = 1
# >alternative hypothesis: significant effects

y,

lmtest::coxtest(dfMo[[1,2]], dfMo[[2,2]])
# > Cox test
# > 
# > Model 1: math ~ size + year
# > Model 2: math ~ size + year
# >                 Estimate Std. Error    z value Pr(>|z|)    
# > fitted(M1) ~ M2     0.32    1.66695     0.1898   0.8494    
# > fitted(M2) ~ M1 -1222.87    0.13616 -8981.1963   <2e-16 ***
# > ---
# > Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# > Warning messages:
# > 1: In lmtest::coxtest(dfMo[[1, 2]], dfMo[[2, 2]]) :
# >   models fitted on different subsets
# > 2: In lmtest::coxtest(dfMo[[1, 2]], dfMo[[2, 2]]) :
# >   different dependent variables specified

En segundo lugar, estoy interesado en comparar los coeficientes de regresión entre dos grupos. Digamos, ¿la estimación year1.5de 3.04 es significativamente diferente de 3.17? Cf. punto 2 arriba.

Pregunte si alguno de los puntos anteriores no está claro y me complacerá explicarlo. ¡Cualquier ayuda será apreciada!

Me doy cuenta de que esta pregunta es un poco de programación, pero inicialmente la publiqué en SO. Sin embargo, DWin tuvo la amabilidad de señalar que la pregunta pertenecía a CrossValidated y la migró aquí.

Eric Fail
fuente
@DWin, gracias. Lo publiqué en SO ya que anteriormente obtuve algunas respuestas realmente buenas con respecto a este tipo de modelos y el plmpaquete, en stackoverflow.com. En el futuro, tendré más cuidado de publicar mis preguntas en el lugar apropiado. Gracias.
Eric Fail
2
No piense que la prueba F funcionaría aquí, ya que sus dos modelos actuales (femenino y masculino) no están anidados. ¿Por qué no incluir run plm con términos de interacción entre variables femeninas y explicativas, por ej plm(math ~ Female * (x1 + x2)). Para probar la primera hipótesis nula, que acaba de ejecutar la prueba F para todos los coeficientes asociados Female:x1, Female:x2. Para probar el segundo nulo, solo necesita probar el parámetro asociado con Female:year1.5.
Semibruin
1
Gracias por tu comentario. Estoy de acuerdo, con respecto a que la prueba F no sea apropiada aquí. Agradezco su sugerencia, pero tengo que implementar esto en un contexto donde la solución de interacción podría no ser factible. Sin embargo, si tiene tiempo, le sugiero que publique su solución como respuesta. Quizás inspire a otros que tienen un problema similar.
Eric Fail
1
Recientemente también me encontré con este problema, pero no pude resolverlo en R. Usé Stata entonces, donde podemos aplicar suestpara ver si dos modelos son significativamente diferentes. Hay una suest()función en un paquete para R, pero dudo que sea la misma. En Stata suestse relaciona con "Estimación aparentemente no relacionada". Tenga en cuenta que eso sureges algo diferente. También estoy interesado en una solución R. Espero que eso ayude de alguna manera.
jay.sf
1
@jaySf, gracias por tu aporte. Tal vez necesitemos migrar esta pregunta nuevamente a stackoverflow.com para descubrir cómo se hace esto en r . No he usado stata en años. ¿Podrías señalar alguna documentación? Gracias.
Eric Fail

Respuestas:

3

El siguiente código implementó la práctica de poner interacción entre Femaleficticio y año. La prueba F en la parte inferior prueba tu nuloβFmimetrounalmi=βMETROunalmi. La estadística t de plmsalida prueba su nuloβFmimetrounalmi:ymiunar=1,5=βMETROunalmi:ymiunar=1,5. En particular, para year=1.5, el valor p es 0.32.

library(plm)  # Use plm
library(car)  # Use F-test in command linearHypothesis
library(tidyverse)
data(egsingle, package = 'mlmRev')
dta <- egsingle %>% mutate(Female = recode(female, .default = 0L, `Female` = 1L))
plm1 <- plm(math ~ Female * (year), data = dta, index = c('childid', 'year', 'schoolid'), model = 'within')

# Output from `summary(plm1)` --- I deleted a few lines to save space.
# Coefficients:
#                 Estimate Std. Error t-value Pr(>|t|)    
# year-1.5          0.8842     0.1008    8.77   <2e-16 ***
# year-0.5          1.8821     0.1007   18.70   <2e-16 ***
# year0.5           2.5626     0.1011   25.36   <2e-16 ***
# year1.5           3.1680     0.1016   31.18   <2e-16 ***
# year2.5           3.9841     0.1022   38.98   <2e-16 ***
# Female:year-1.5  -0.0918     0.1248   -0.74     0.46    
# Female:year-0.5  -0.0773     0.1246   -0.62     0.53    
# Female:year0.5   -0.0517     0.1255   -0.41     0.68    
# Female:year1.5   -0.1265     0.1265   -1.00     0.32    
# Female:year2.5   -0.1465     0.1275   -1.15     0.25    
# ---

xnames <- names(coef(plm1)) # a vector of all independent variables' names in 'plm1'
# Use 'grepl' to construct a vector of logic value that is TRUE if the variable
# name starts with 'Female:' at the beginning. This is generic, to pick up
# every variable that starts with 'year' at the beginning, just write
# 'grepl('^year+', xnames)'.
picked <- grepl('^Female:+', xnames)
linearHypothesis(plm1, xnames[picked])

# Hypothesis:
# Female:year - 1.5 = 0
# Female:year - 0.5 = 0
# Female:year0.5 = 0
# Female:year1.5 = 0
# Female:year2.5 = 0
# 
# Model 1: restricted model
# Model 2: math ~ Female * (year)
# 
#   Res.Df Df Chisq Pr(>Chisq)
# 1   5504                    
# 2   5499  5  6.15       0.29
semibruin
fuente
Muy interesante. Lo intentaré con mis datos de producción. Gracias. Puede publicar la misma respuesta aquí stackoverflow.com/questions/28334298/… y obtener la recompensa allí también.
Eric Fail
Pregunta rápida, ¿crees que es posible reescribir el -c(1:5)bloque de alguna manera que haga que el código sea más genérico? Tengo vectores de tamaño cambiante entrando y saliendo y una respuesta más genérica también podría beneficiar a otros.
Eric Fail
@EricFail Reemplacé -c(1:5)con expresión regular. Es más genérico ahora. En general, le gustaría usar greplpara unir patrones en presencia de muchas variables.
Semibruin