Waldtest in R to get adjust F statistic with plm and result shown with stargazer?

Question

I am working with an unbalanced short panel. Raw data: bankFull.xlsx

What I actually want is only get the regression results with two side fixed effects and robust S.E reported, which is very easy in Stata. I followed online tutorial but ran into some problem always with

# Adjust F statistic 
wald_results <- waldtest(FE1, vcov = cov1)
Error in model.matrix.pFormula(formula, data, rhs = 1, model = model,  : 
  NA in the individual index variable

no matter how I adjusted the data! It almost drives me crazy.

here is my code:

bankFull <- openxlsx::read.xlsx("bankFull.xlsx",1)

attach(bankFull)
library(plm)

FE1 = plm(  RoA ~
              log(1+degreeNW)+
              ln_assets+
              log(no_of_board_members/staffNo)+
              log(no_of_branch_covered_city)+ 
              log(operation_year)+
              `RoA-1`+
              log(staffNo),
            data = bankFull, index = c("name","year"),  
            effect="twoways",na.action = na.omit,
            model= "within")

# robust S.E.-----------
library(sandwich)
library(lmtest)   # waldtest; see also coeftest.
library(stargazer)

# Adjust standard errors
cov1         <- vcovHC(FE1, type = "HC1")
robust_se    <- sqrt(diag(cov1)) 

# Adjust F statistic 
wald_results <- waldtest(FE1, vcov = cov1)

# show results. how can I get the F value?
stargazer(FE1, FE1, type = "text",
          se        = list(NULL, robust_se),
          omit.stat = "f")

Secondly, as the code shown, I use stargazer to demonstrate the results. I also need the adjusted F value to be shown in the table. Is there any option in the package that I can use?

Answer 1

Edit: update information according to CRAN release 1.6-4 of plm

Use the CRAN version 1.6-4 of plm supports robust F tests for your model by function pwaldtest (this function was called Ftest in the development version but renamed prior to CRAN release).

Example:

data("Grunfeld", package = "plm")
mod_fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within")
plm::pwaldtest(mod_fe, test = "F")


# with robust vcov
plm::pwaldtest(mod_fe, test = "F", vcov = vcovHC(mod_fe))
plm::pwaldtest(mod_fe, test = "F", vcov = function(x) vcovHC(x, type = "HC3"))
summary(mod_fe, vcov = vcovHC)

To feed the robust values (robust standard errors, t- and p-values, F-value associated p-value) use arguments se, t, p and for the F test simply add.lines of the stargazer command (and omit the F statistic generated by stargazer by default). Here is a full example for what you want: http://jakeruss.com/cheatsheets/stargazer.html (section "Robust standard errors (replicating Stata’s robust option)").