I am performing an OLS regression with time and industry fixed effects and I want to compare a potential selection problem with the Heckman approach in R.

In R, I have the following code for my OLS regression:

# Load necessary libraries
library(plm)
library(sampleSelection)

# OLS Model with fixed effects - with complete control variables
ols.model <- plm(VALUE ~ Score + SIZE + LEVERAGE + ROA + AGE + GROWTH, 
                         data = dataset,
                         model = "within",
                         index = c("Year", "NAICS"))
summary(ols.model) 

# Run Heckman selection model
selection_eg <- as.formula(SELECTION_VARIABLE ~ SIZE + INDUSTRY + AVAILABILITY)

outcome_eq <- as.formula(VALUE ~ Score + SIZE + LEVERAGE + ROA + AGE + GROWTH)

heckman_model <- selection(selection_eq, 
                           outcome_eq, 
                           data = dataset, 
                           method = "2step")

summary(heckman_model)

My question: How can I include the time and industry fixed effects like I did in the OLS regression index = c(“Year”, “NAICS”)?

I tried including the index = c(“Year”, “NAICS”) in the Heckman function, but this is not working due to the structure of the formula.