Treatment effects Coefplot: Please help

[Minda Bel]

Hello Statalist, please help!
I am trying to replicate the coefplot figure attached below
I have created a sample data for you and shared my code. please help.

*generate sample data with 100obs first.
clear
set obs 100

* Generate binary treatment variable (random 0 or 1)
gen treatment = round(runiform())

* Generate control variables (random normal distribution)
gen control1 = rnormal()
gen control2 = rnormal()

* Generate gender variable (50% Male, 50% Female)
gen female = round(runiform())

* Generate outcome variables based on treatment & controls with random noise
gen overall_index = 0.3*treatment + 0.2*female + rnormal()
gen intent_stay = 0.25*treatment - 0.1*female + rnormal()
gen job_search = -0.2*treatment + 0.15*female + rnormal()
gen hours_worked = 0.1*treatment - 0.05*female + rnormal()
gen employment = 0.35*treatment + 0.05*female + rnormal()
gen labor_income = 0.4*treatment - 0.1*female + rnormal()
gen food_security = 0.2*treatment + 0.1*female + rnormal()
gen consumption = 0.15*treatment + 0.05*female + rnormal()
gen life_satisfaction = 0.3*treatment - 0.05*female + rnormal()
gen whodas_score = -0.25*treatment + 0.1*female + rnormal()


* Define the list of outcome variables
local outcomes overall_index intent_stay job_search hours_worked employment ///
labor_income food_security consumption life_satisfaction whodas_score

* Create a temporary dataset to store results
tempname memhold
postfile `memhold' str20 outcome str10 group coef se using myresults.dta, replace

* Loop through each outcome variable and estimate separate regressions for males and females
foreach outcome of local outcomes {
* Regression for males (female == 0)
quietly regress `outcome' treatment $con if female == 0
post `memhold' ("`outcome'") ("Male") (_b[treatment]) (_se[treatment])

* Regression for females (female == 1)
quietly regress `outcome' treatment $con if female == 1
post `memhold' ("`outcome'") ("Female") (_b[treatment]) (_se[treatment])
}

postclose `memhold'
* Load the dataset
use myresults.dta, clear

* Check for duplicate entries
duplicates report outcome group
duplicates list outcome group

* If duplicates exist, remove them before reshaping
duplicates drop outcome group, force

* Reshape the dataset for coefplot
reshape wide coef se, i(outcome) j(group) string

* Generate plot using coefplot
coefplot (coefMale, label("Male") color(purple) msymbol(diamond)) ///
(coefFemale, label("Female") color(orange) msymbol(diamond)), ///
vertical ///
yline(0, lcolor(black)) ///
legend(order(1 "Male" 2 "Female")) ///
xtitle("Effect Size") xlabel(-0.4(0.2)0.8) ///
coeflabels(overall_index = "Overall Index" ///
intent_stay = "Intent to Stay" ///
job_search = "Job Search" ///
hours_worked = "Hours Worked" ///
employment = "Employment" ///
labor_income = "Labor Income" ///
food_security = "Food Security" ///
consumption = "Consumption" ///
life_satisfaction = "Life Satisfaction" ///
whodas_score = "WHODAS Score") ///
bycoefs ///
grid

[Andrew Musau]

coefplot is from SSC (FAQ Advice #12). Thanks for the reproducible example. Your results are in a dataset, so you can partition them into several matrices and plot.

*generate sample data with 100obs first.
clear
set obs 100

* Generate binary treatment variable (random 0 or 1)
gen treatment = round(runiform())

* Generate control variables (random normal distribution)
gen control1 = rnormal()
gen control2 = rnormal()

* Generate gender variable (50% Male, 50% Female)
gen female = round(runiform())

* Generate outcome variables based on treatment & controls with random noise
gen overall_index = 0.3*treatment + 0.2*female + rnormal()
gen intent_stay = 0.25*treatment - 0.1*female + rnormal()
gen job_search = -0.2*treatment + 0.15*female + rnormal()
gen hours_worked = 0.1*treatment - 0.05*female + rnormal()
gen employment = 0.35*treatment + 0.05*female + rnormal()
gen labor_income = 0.4*treatment - 0.1*female + rnormal()
gen food_security = 0.2*treatment + 0.1*female + rnormal()
gen consumption = 0.15*treatment + 0.05*female + rnormal()
gen life_satisfaction = 0.3*treatment - 0.05*female + rnormal()
gen whodas_score = -0.25*treatment + 0.1*female + rnormal()


* Define the list of outcome variables
local outcomes overall_index intent_stay job_search hours_worked employment ///
labor_income food_security consumption life_satisfaction whodas_score

* Create a temporary dataset to store results
tempname memhold
postfile `memhold' str20 outcome str10 group coef se using myresults.dta, replace

* Loop through each outcome variable and estimate separate regressions for males and females
foreach outcome of local outcomes {
* Regression for males (female == 0)
quietly regress `outcome' treatment $con if female == 0
post `memhold' ("`outcome'") ("Male") (_b[treatment]) (_se[treatment])

* Regression for females (female == 1)
quietly regress `outcome' treatment $con if female == 1
post `memhold' ("`outcome'") ("Female") (_b[treatment]) (_se[treatment])
}

postclose `memhold'
* Load the dataset
use myresults.dta, clear

* Check for duplicate entries
duplicates report outcome group
duplicates list outcome group

* If duplicates exist, remove them before reshaping
duplicates drop outcome group, force

* Reshape the dataset for coefplot
reshape wide coef se, i(outcome) j(group) string

mkmat coefFemale seFemale if _n==1, mat(overallF) rownames(outcome)
mkmat coefMale seMale if _n==1, mat(overallM) rownames(outcome)

mkmat coefFemale seFemale if inrange(_n,2, 5), mat(targetF) rownames(outcome)
mkmat coefFemale seMale if inrange(_n,2, 5), mat(targetM) rownames(outcome)

mkmat coefFemale seFemale if inrange(_n,6, 10), mat(downF) rownames(outcome)
mkmat coefFemale seMale if inrange(_n,6, 10), mat(downM) rownames(outcome)


* Generate plot using coefplot
coefplot (matrix(overallF[,1]), mc(red) ciopts(lc(red)) se(overallF[,2]) aseq("Mental Experiencing")) ///
(matrix(overallM[,1]), mc(blue) ciopts(lc(blue) mc(blue)) se(overallM[,2]) aseq("Mental Experiencing")) ///
(matrix(targetF[,1]), mc(red) ciopts(lc(red)) se(targetF[,2]) aseq("Targeted Behaviors")) ///
(matrix(targetM[,1]),  mc(blue) ciopts(lc(blue) mc(blue)) se(targetM[,2]) aseq("Targeted Behaviors")) ///
(matrix(downF[,1]), mc(red) ciopts(lc(red)) se(downF[,2]) aseq("Downstream Outcomes")) ///
(matrix(downM[,1]), mc(blue) ciopts(lc(blue) mc(blue)) se(downM[,2]) aseq("Downstream Outcomes")), eqstrict ///
eqlabel(`""Mental     " "Experiencing""' `""Targeted " "Behaviors""' `""Downstream" "Outcomes  ""') ///
legend(order(2 "Female" 4 "Male")) ///
xtitle("Effect Size") xlabel(-0.4(0.2)0.8) ///
coeflabels(overall_index = "Overall Index" ///
intent_stay = "Intent to Stay" /// 
job_search = "Job Search" ///
hours_worked = "Hours Worked" ///
employment = "Employment" ///
labor_income = "Labor Income" ///
food_security = "Food Security" ///
consumption = "Consumption" ///
life_satisfaction = "Life Satisfaction" ///
whodas_score = "WHODAS Score")