Help with Implementing Hayes' Mediation Analysis in STATA

[gi peters]

Hi,

Im trying to implement mediation analysis using Dr. Hayes' approach in STATA and has been having some issues with how long it takes. One of my codes has been running for 4 days and has only managed to do 50 replications at this point. I do my analysis in a data center so I couldnt change anything when it comes to computer hardware. Could anyone be kind to comment on my code and see if I could do something about it to expedite the process?


This is my code for setting the survey parameters and determining the a, b and c' paths and the legend for these coefficients

svyset [pweight=surveyweight], bsrweight(weight1- weight500) vce(bootstrap) mse
local controls "cc_agesq ib1.cc_sex ib2.cc_educat4 ib1.cc_race5 ib1.cc_marstat4 c.cc_health ib0.cc_chronic3 ib1.cc_language ib1.wave_codes"
svy:gsem (boyd_scores `controls'->cc_combelong)(boyd_scores `controls'->cc_personalincome)(cc_combelong cc_personalincome boyd_scores `controls'->cc_mhealth) if flag1==1 & immigrant3==0 
gsem, coeflegend

Now, I want to determine indirect effects through 1)combelong and 2)income; total indirect effects; and total effects. I saw in the STATA manual that svy: works with nlcom so I tried but it did not work. This is my code that did not work

svy: nlcom_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong]
svy: nlcom _b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome]
svy: nlcom (_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong])+(_b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome])    
svy: nlcom _b[cc_mhealth:boyd_scores]+((_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong])+(_b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome]))

I searched online and found a code here in STATAlist but the post is from 2020 and when I posted my question there, the original author of the question did not respond. Maybe no longer active in STATAlist. This code is working but it's been 4 days and it only managed to do 50 replications. Could this be improved or is there another code that I could use that could provide the same output?

capture program drop boot1
program boot1, rclass
local controls "cc_agesq ib1.cc_sex ib2.cc_educat4 ib1.cc_race5 ib1.cc_marstat4 c.cc_health ib0.cc_chronic3 ib1.cc_language ib1.wave_codes"
svy: gsem (boyd_scores `controls'->cc_combelong)(boyd_scores `controls'->cc_personalincome)(boyd_scores cc_combelong cc_personalincome `controls'->cc_mhealth) if flag1==1 & immigrant3==0 
return scalar indirectcom=(_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong])
return scalar indirectinc=(_b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome])
return scalar indirecttot=(((_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong])+(_b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome])))    
return scalar effecttot=(_b[cc_mhealth:boyd_scores]+((_b[cc_combelong:boyd_scores]*_b[cc_mhealth:cc_combelong])+(_b[cc_personalincome:boyd_scores]*_b[cc_mhealth:cc_personalincome])))
end 
bootstrap r(indirectcom) r(indirectinc) r(indirecttot)r(effecttot), seed(123) reps(500): boot1 estat boot, bc percentile

Thanks in advance!