Probit classification models with country fixed effects (panel data) in R

I want to do a probit regression for a panel data with country fixed effects and Clustered (by country) standard errors.

The dependent variable is y and the independent the x. The panel data contains the year and the country.

I have tried to solve this problem as follows:

     library(pglm)
     model <- plm( y ~ x1 + x2 + factor(country), data = c1, family = binomial(link = "probit"))
        summary(model)

But I didn't get the values I expect. As a reference I have a STATA code, where in contrast to my idea, there were created country dummy variables. Maybe that would be also a method.

Thank you for your help!