How can I add a multidimensional panel to estimate TFP (year, id, region) in R?

Question

I'm not sure about adding fixed effects for the variables year,id,region (using the estprod library) as in this paper using LP stimator.

In the STATA forum ( #18 #19; prodest library) it seems that to add fixed effects I have to create dummies Control parameter.

qui tab year, gen(dy) // generate one dummy per year

prodest log_y, free(log_lab1 log_lab2) state(log_k) proxy(log_investment) va met(op) poly(4) reps(40) id(id) t(year) control(dy*)

To make a third panel PLM I found in STATA as in R

Code

library(dplyr)
library(estprod)

data = read.csv("data_tfp.csv")

data = data %>% 
  mutate(v = 1,
         rgn = region,
         group = group_indices(., id, region),
         grp = group) %>% 
  spread(rgn, v, fill = 0) %>% #generate one dummy per region
  arrange(id, year)

md0 = estprod::levinsohn_petrin(data,
                                formula = y ~ l | k | m | a+b+c+d,
                                id = "id",
                                time = "year",
                                reps = 10,
                                gross = T)
summary(md0)

My data

here