Comparing mean values by group between several variables

I am trying to reproduce a graph from Stata in R. I have several variables and want to display their mean in each treatment group of which there are two. The Stata graph is as follows:

This coefficient plot is not actually a plot of coefficients, but of the mean values by each treatment for each separate variable. The df basically looks something like.

workable data

I figured out geom_pointrange() is probably what you are looking for:

library("ggplot2")
set.seed(111018)
interval1 <- -qnorm((1-0.9)/2)  

means_treatment_1 <- rnorm(2)
se_treatment_1 <- rnorm(2)

df_treatment_1 <- data.frame("Mean" = means_treatment_1,
                         "lower" = means_treatment_1 - se_treatment_1*interval1,
                         "upper" = means_treatment_1 + se_treatment_1*interval1,
                         "Variable" = c("medicare_spending_dummy", 
                                        "job_training_dummy"),
                         "Treatment" = "a")


means_treatment_2 <- rnorm(2)
se_treatment_2 <- rnorm(2)

df_treatment_2 <- data.frame("Mean" = means_treatment_2,
                         "lower" = means_treatment_2 - se_treatment_2*interval1,
                         "upper" = means_treatment_2 + se_treatment_2*interval1,
                         "Variable" = c("medicare_spending_dummy", 
                                        "job_training_dummy"),
                         "Treatment" = "b")



df_tot<-rbind(df_treatment_1, df_treatment_2)



# Plot


ggplot(df_tot, aes(colour = Treatment)) +
geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
geom_pointrange(aes(x = Variable, y = Mean, ymin = lower, ymax = upper ),lwd = 1, position = position_dodge(width = 1/2)) +

coord_flip() + 
theme_bw()