Can I add 3 fixed effect (firm, industry, year) in one Fixed effect panel model in R? I am using the formula below to estimate the fixed effect panel model in R:
model_fe <- plm(formula = y~ x, data, 
                index = c("firm ID", "year"), 
                model = "within", effect = "individual")

However, I need to include 3 fixed effects simultaneously in the model (firm, industry, and year), but not sure how to do that in R. Can I include the 'industry' term in the 'index' function? If yes, what will be the effect? two ways? Please suggest if there is also some other way to simultaneously accommodate all the three.
I think I can also use dummy variables by using the function factors(), but that will add additional variables to the model which will in turn decrease the degree of freedom and inflate the r-squared.

It seems you're dealing with multidimensional panel data. Maybe this link that can help you.
However, if the time variable is $year$, the group variable should be $firm$ and the extra-dimension will be $industry$. How many industries/year/industries do you have? The size of each dimension is critical.
The problem with multiple effects models (or more) is that practically you're doing a model for each exact dimension (industry-firm-year) with few degrees of freedom each. Wooldridge has some literature about that.
Perhaps, you need to review why is that you need the multiple effects. To view results on that level or for statistical treatment purposes. Too many factors could create an overfitted model.
