Extract fixed effects and match with values of time invariant variables, for testing time invariant effects.

fixef_of_time_invariant_vars.R

## Function to extract fixed effects from plm or fixest model, and match with values of time invariant variables.
## Arguments:
##  model = a plm or fixest fixed effects model object
##  factors = either a data.frame/pdata.frame/data.table object containing time-invariant variables and a panel index variable,
##            or in the case that the model arg is a plm object, factors can be a vector of names of time-invariant variables
##            that will be extracted from the model object (along with the panel index).
##  idvar = the name of the panel index variable; by default this will be set to the name of the model index if model is
##          a plm object, or the name of the factors index if factors is a pdata.frame index. Otherwise if model is a fixest
##          object then idvar must specify a variable in factors (a data.frame object).
## For example of usage see code after this function.
fixef_of_time_invariant_vars <- function(model, factors, idvar = NA) {
    require(data.table)
    fes <- fixef(model)
    if("plm" %in% class(model)) {
        modeldata <- model$model
        if(is.na(idvar)) idvar <- names(index(model))[1]
        fes.dt <- data.table(fe = fes)
        fes.dt[,(idvar) := names(fes)]
        setkeyv(fes.dt, idvar)
        idxs <- index(model$model)[1]
        if(is.character(factors)) {
            modeldata.dt <- as.data.table(model$model)
            modeldata.dt[,(idvar) := idxs]
            setkeyv(modeldata.dt, idvar)
            dt1 <- modeldata.dt[,.SD,.SDcols=c(factors,idvar)]
        } else if("data.frame" %in% class(factors)) {
            dt1 <- as.data.table(factors,key = idvar)
        } else stop("factors arg must be a list of factor names, or a data.frame")
    } else if("fixest" %in% class(model)) {
        stopifnot("data.frame" %in% class(factors))        
        if("pdata.frame" %in% class(factors)) {
            if(is.na(idvar)) {
                idvar <- names(index(factors))[1] 
            }
        } else stopifnot(!is.na(idvar))
        idxs <- names(fes[[1]])
        fes.dt <- data.table(fe = fes[[1]])
        fes.dt[,(idvar) := idxs]
        setkeyv(fes.dt, idvar)
        ## TODO: should I be storing the indices for the other fixed effects?
        if(length(fes)>1) {
            for(i in 2:length(fes)) {
                dt1 <- data.table(fe = fes[[i]])
                dt1[,(idvar) := idxs]
                setkeyv(dt1, idvar)
                fes.dt <- dt1[fes.dt]
            }
        }
        if("data.frame" %in% class(factors)) {
            dt1 <- as.data.table(factors,key = idvar)
        } else stop("factors arg must be a data.frame object")
    }
    dt1 <- dt1[,lapply(.SD,first),by=idvar]
    fes.dt <- dt1[fes.dt]
    return(fes.dt)
}

This function is useful when you are doing panel data modelling and need to test a time invariant effect, but the Hausmann test is significant which suggests you should use a fixed effects model (which doesn't allow time invariant regressors).
It works with plm and fixest models in R.
To test the significant of the time invariant variables, first extract the fixed effects using the function defined in this gist, then estimate a regression model with the fixed effects as dependent variable, and the time invariant variable(s) as independent variable(s), e.g:

Extract fixed effects and corresponding values of CS_C1:
tivdata <- fixef_of_time_invariant_vars(plm_model,"CS_C1")
Examine effect of CS_C1 on fixed effects of plm_model:
lm(fe~I(as.factor(CS_C1)),data = tivdata)