8
$\begingroup$

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    
#> Adj. R^2      0.70         0.72    
#> 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.

$\endgroup$
  • $\begingroup$ @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. $\endgroup$ – Eric Fail Jun 17 '18 at 8:10
  • 2
    $\begingroup$ 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. $\endgroup$ – semibruin Jun 18 '18 at 22:50
  • 1
    $\begingroup$ 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. $\endgroup$ – Eric Fail Jun 19 '18 at 7:03
  • 1
    $\begingroup$ 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. $\endgroup$ – jay.sf Jun 19 '18 at 13:50
  • 1
    $\begingroup$ @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. $\endgroup$ – Eric Fail Jun 19 '18 at 14:19
3
+50
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Eric Fail Jun 19 '18 at 20:29
  • $\begingroup$ 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. $\endgroup$ – Eric Fail Jun 20 '18 at 4:23
  • $\begingroup$ @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. $\endgroup$ – semibruin Jun 20 '18 at 4:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.