comparing groups in repeated measures FE models, with a nested error component, estimated using plm

I have estimated some repeated measures Fixed Effects models, with a nested error component, based on on grouping variables, i.e. non-nested models, using . I am now interested to

1. test if the full models are significantly different, i.e. $$H_o: \beta_{Female} = \beta_{Male}$$ where $\beta_{Female}$ is the full model for Females and $\beta_{Male}$ is the full model for Males and
2. subsequently test selected regression coefficients between two groups, i.e. $$H_o: \beta_{Female == year1.5} = \beta_{Male == year1.5}$$ where $\beta_{Female == year1.5}$ is the regression coefficient for females at year1.5, and $\beta_{Male == year1.5}$ is the regression coefficient for males at year1.5.

I will illustrate the situation using the below working example,

First, some packages needed,

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


Second, some data preparation,

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


Third, I estimate a set of models for each gender in data

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


Forth, lets look at the two estimated models,

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
#> Num. obs.  3545         3685
#> ===================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05    #>


Now, I want to test if these two (linear OLS) models are significantly different, cf. point1 above. I looked around SO and the internet and some suggest that I need to use plm::pFtest(), also suggested here, which I have tried, but I'm not convinced. I would have imagined some test for non-nested models, possible Cox test, lmtest::coxtest, but I am not sure at all. If someone here could possibly help me.

I tried,

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


and,

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


Second, I am interested to compare regression coefficients between two groups. Say, is the estimate for year1.5 of 3.04 significantly different from 3.17? Cf. point 2 above.

Please ask if any of the above is not clear and I will be happy to elaborate. Any help will be greatly appreciated!

I realize this question is a bit programming like, but I initially posted it in SO. However, DWin was kind enough to point out that the question belonged in CrossValidated and migrated it here.

• @DWin, Thanks. I posted it in SO as I've previously gotten some really good answers regarding this types of models, and the plm package, at stackoverflow.com. I will take more care in the future to post my questions in the appropriate place. Thanks. – Eric Fail Jun 17 '18 at 8:10
• Don't think F test would work here, since your current two models (female and male) are not nested. Why not include run plm with interaction terms between female and explanatory variables, e.g. plm(math ~ Female * (x1 + x2)). To test the first null hypothesis, you just run F test for all coefficients associated with Female:x1, Female:x2. To test the second null, you just need t test the parameter associated with Female:year1.5. – semibruin Jun 18 '18 at 22:50
• Thank you for your comment. I agree, with regard to the F-test not being appropriate here. I appreciate your suggestion, but I have to implement this in a context where the interaction solution might not be feasible. However, if you have the time, I suggest you post your solution as an answer. Maybe it will inspire others who has a similar issue. – Eric Fail Jun 19 '18 at 7:03
• I recently came about this issue, too, but couldn't solve it in R. I used Stata then, where we can apply suest to see if two models are significantly different. There's a suest() function around in a package for R but I doubt that it is the same. In Stata suest is related to "Seemingly unrelated estimation". Note, that sureg is somewhat different. I am also interested in an R solution. Hope that would help somehow. – jay.sf Jun 19 '18 at 13:50
• @jaySf, thank you for your input. Maybe we need to migrate this question back to stackoverflow.com to figure out how this is done in r. I haven't used stata in years. Could you possibly point to some documentation? Thanks. – Eric Fail Jun 19 '18 at 14:19

The following code implemented the practice of putting interaction between Female dummy and year. The F test at the bottom test your null $\beta_{Female} = \beta_{Male}$. The t-statistic from plm output tests your null $\beta_{Female:year=1.5}=\beta_{Male:year=1.5}$. In particular, for year=1.5, the p-value is 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

• Very interesting. I'll try it on my production data. Thanks. You could post the same answer here stackoverflow.com/questions/28334298/… and get the bounty there too. – Eric Fail Jun 19 '18 at 20:29
• Quick question, do you think it is possible to rewrite the -c(1:5) block in some way that make the code more generic? I have shifting size vectors going in and out and a more generic answer would possible also benefit others. – Eric Fail Jun 20 '18 at 4:23
• @EricFail I replaced -c(1:5) with regular expression. It is more generic now. In general, you would like to use grepl to match patterns in the presence of a lot of variables. – semibruin Jun 20 '18 at 4:47