Двусторонняя регрессия панели с фиксированными эффектами по группам

Я искал много вопросов о том, как я могу успешно это сделать, но безрезультатно. Самое близкое, что я нашел, это: Регрессия панели по группам , но даже это не полностью охватывает объем того, что я хочу.

У меня есть набор данных, содержащий финансовые данные о нескольких фирмах в разных отраслях за 11 лет. Теперь я успешно выполнил двухстороннюю регрессию с фиксированными эффектами по всему набору данных, используя:

plm(ROS ~ ELI + ELI*Pandemic + Growth + Assets,
               data = df,
               index = c("Year", "Firm"),
               model = "within",
               effect = "twoways")

И у меня довольно хорошие результаты, за исключением того, что они не на 100% соответствуют моим ожиданиям. Теперь я хотел бы провести ту же самую регрессию, но сгруппированную по отраслям (NAICS), в которых работают фирмы. Таким образом, я могу глубоко изучить, имеют ли разные отрасли разные результаты.

Образец моего набора данных можно найти ниже, хотя я думаю, что для успешного выполнения этого требуется больше данных. Если это так, я бы с радостью отправил весь свой набор данных. По сути, у меня есть 46 отраслей (NAICS), в каждой из которых несколько фирм, и данные по каждой фирме за 11 лет.

df <- structure(list(NAICS = c(315, 3151, 3152, 3153, 3154, 3155), 
Year = c(2010, 2011, 2012, 2013, 2014, 2015), 
Firm = c("A", "B", "C", "D", "E", "F"), 
Assets = c(399.5, 674.42, 1289.57, 2216.97, 2691.89, 2566.8), 
Inventories = c(117.17, 187.41, 266.89, 426.94, 519.91, 546.8), 
Sales = c(803.34, 1302.25, 2181.73, 3310.84, 4371.47, 4712.1), 
Growth = c(0, -0.62, -0.68, -0.52, -0.32, -0.08), 
ROS = c(0.17, 0.19, 0.29, 0.3, 0.29, 0.25), 
ELI = c(0.856264428748774, 0.723379402777553, 0.958341156943977, 0.680567730897854, 0.790480861209701, 0.827279134948296), 
Pandemic = c(0, 0, 0, 0, 0, 0)), 
row.names = c(NA, -6L), 
class = c("tbl_df", "tbl", "data.frame"))

Я надеюсь, что это достаточно ясно.