I have to replicate a Paper with a dprobit regression from STATA to R.
Unfortunately, my estimates in R differ significantly from the STATA ones (i.e. 0.044 in STATA and 0.125 in R for productivity_indicator3)
Does anyone know what to bear in mind when replicating a probit regression in R?
Attached you can find one example from the STATA code and my corresponding Try with R:

STATA:

dprobit piece_rate productivity_indicator3 risk_attitude self_assessment trust coeff_trans         female, robust cluster(session_id)
outreg2 using "table2.xls", replace tex excel label se br bdec(3) rdec(3) ctitle(Piece rate) addstat(Pseudo R-squared, `e(r2_p)')

R:

dat = data%>%  
mutate(pr = ifelse(piece_rate == "Piece Rate", 1,0)) %>%  


regression_pr = glm(pr ~ productivity_indicator3 + risk_attitude + self_assessment + trust   + coeff_transfer1_nocon + female, data = dat, family = binomial("probit"))

cluster_se = vcovCL(regression_pr, cluster = ~ session_id)
coeftest(regression_pr, vcov = cluster_se)

fyi: this is about an experiment where there are three treatments. The code above covers just one treatment as an example. The initial data provides the variable piece_rate with the two corresponding values c(“Piece Rate”, “Fix Payment”). I have created a new variable from this (“pr”) which is for every treatment == Piece Rate = 1 and otherwise 0, to prepare the probit regression.
There are also no missings in the data frame.

Would be very happy about any idea what I did wrong or what to keep in mind when doing this. Thanks in advance! :))

I tried the following, but the estimates differ significantly from those in STATA

dat = data%>%  
mutate(pr = ifelse(piece_rate == "Piece Rate", 1,0)) %>%  


regression_pr = glm(pr ~ productivity_indicator3 + risk_attitude + self_assessment + trust + coeff_transfer1_nocon + female, data = dat, family = binomial("probit"))

cluster_se = vcovCL(regression_pr, cluster = ~ session_id)

coeftest(regression_pr, vcov = cluster_se)