comparando grupos em medidas repetidas modelos FE, com um componente de erro aninhado, estimado usando plm

8

Estimei algumas medidas repetidas modelos de efeitos fixos, com um componente de erro aninhado, com base em variáveis ​​de agrupamento, ou seja, modelos não aninhados, usando . Agora estou interessado em

  1. teste se os modelos completos são significativamente diferentes, ou seja, onde é o modelo completo para e é o modelo completo para e β F e m a l e β M a l e
    Ho:βFemumaeue=βMumaeue
    βFemumaeueFemalesβMumaeueMales
  2. posteriormente teste os coeficientes de regressão selecionados entre dois grupos, ou seja, que é o coeficiente de regressão para mulheres at , e é o coeficiente de regressão para homens em . β F e m um l e = = Y e um r 1,5 β M um l e = = y e a r 1,5
    Ho:βFemumaeue==yeumar1.5=βMumaeue==yeumar1.5
    βFemumaeue==yeumar1.5year1.5βMumaeue==yeumar1.5year1.5

Ilustrarei a situação usando o exemplo de trabalho abaixo,

Primeiro, alguns pacotes necessários,

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

Segundo, alguma preparação de dados,

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

Terceiro, eu estimo um conjunto de modelos para cada gênero nos dados

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") )

Quarto, vamos olhar para os dois 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    #> 

Agora, quero testar se esses dois modelos (OLS lineares) são significativamente diferentes, cf. ponto 1 acima. Eu olhei em volta da SO e da Internet e alguns sugerem que eu preciso usar plm::pFtest(), também sugerido aqui , o que eu tentei, mas não estou convencido. Eu teria imaginado algum teste para modelos não aninhados, possível teste de Cox lmtest::coxtest, mas não tenho certeza. Se alguém aqui pudesse me ajudar.

Eu tentei,

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

e,

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

Segundo, estou interessado em comparar coeficientes de regressão entre dois grupos. Digamos, a estimativa year1.5de 3,04 é significativamente diferente de 3,17? Cf. ponto 2 acima.

Por favor, pergunte se alguma das opções acima não está clara e terei prazer em elaborar. Qualquer ajuda será muito apreciada!

Sei que essa pergunta é um pouco de programação, mas inicialmente a publiquei no SO. No entanto, o DWin teve a gentileza de apontar que a pergunta pertencia ao CrossValidated e a migrou para aqui.

Eric Fail
fonte
@ DWin, obrigado. Eu publiquei no SO, pois já obtive respostas realmente boas sobre esse tipo de modelo e o plmpacote no stackoverflow.com. Terei mais cuidado no futuro para postar minhas perguntas no local apropriado. Obrigado.
Eric falha
2
Não pense que o teste F funcionaria aqui, pois seus dois modelos atuais (feminino e masculino) não estão aninhados. Por que não incluir run plm com termos de interação entre variáveis ​​femininas e explicativas, por exemplo plm(math ~ Female * (x1 + x2)). Para testar a primeira hipótese nula, basta executar o teste F para todos os coeficientes associados Female:x1, Female:x2. Para testar o segundo nulo, você só precisa testar o parâmetro associado Female:year1.5.
Semibruin
1
Obrigado pelo seu comentário. Concordo que o teste F não é apropriado aqui. Agradeço sua sugestão, mas preciso implementá-la em um contexto em que a solução de interação possa não ser viável. No entanto, se você tiver tempo, sugiro que você publique sua solução como resposta. Talvez isso inspire outras pessoas que têm um problema semelhante.
Eric Falha
1
Recentemente, também resolvi esse problema, mas não consegui resolvê-lo na R. Eu usei o Stata então, onde podemos aplicar suestpara ver se dois modelos são significativamente diferentes. Existe uma suest()função em um pacote para R, mas duvido que seja o mesmo. Em Stata, suestestá relacionado à "estimativa aparentemente não relacionada". Note que isso suregé um pouco diferente. Também estou interessado em uma solução R. Espero que ajude de alguma forma.
Jay.sf
1
@jaySf, obrigado pela sua contribuição. Talvez seja necessário migrar essa pergunta de volta para stackoverflow.com para descobrir como isso é feito em r . Eu não uso o stata há anos. Você poderia apontar para alguma documentação? Obrigado.
Eric Falha

Respostas:

3

FemaleβFemumaeue=βMumaeueplmβFemumaeue:yeumar=1.5=βMumaeue:yeumar=1.5year=1.5, o valor p é 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
semibruína
fonte
Muito interessante. Vou tentar nos meus dados de produção. Obrigado. Você pode postar a mesma resposta aqui stackoverflow.com/questions/28334298/… e obter a recompensa lá também.
Eric Falha
Pergunta rápida, você acha que é possível reescrever o -c(1:5)bloco de alguma maneira que torne o código mais genérico? Eu tenho vetores de mudança de tamanho entrando e saindo e uma resposta mais genérica possível também beneficiaria outros.
Eric Falha
@EricFail Eu substituí -c(1:5)pela expressão regular. Agora é mais genérico. Em geral, você gostaria de usar greplpara corresponder padrões na presença de muitas variáveis.
Semibruin