BOOTSTRAP: How to bootstrap iterating many procedures

[saramanzella]

Hello everyone!

I'm using SAS Studio.

I'm writing my master thesis and I'm doing some mediation analysis with time-dependent exposure and mediator.

The code I'm using is the one at the end of the post and when running it, it produces the estimates of direct, indirect and total effects saved in a dataset as follows:

      DIRECT            INDIRECT           TOTAL

0.200344234

-0.002250784

0.1980934497

At this point what I have left to calculate are BOOTSTRAP CONFIDENCE INTERVALS fore those effects estimates.

The only problem is that I've never done bootstrap on SAS iterating such a long code and saving at the end in the same dataset all the estimates and use that as a distribution to get the 5th and 95th percentile as my confidence boundaries.

I was wondering if anyone had any idea on how I can do that or somewhere I can find this information.

Many thanks,

Sara

* TIME-DEPENDENCE ;
* MSM Y ;
* Numerator Weights for the Exposure and Mediator in the Outcome MSM ;
proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w2 = ricchezza_factor_w1 fisico_w1;
	output out = finale.all_stdz student = rn2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w2 = ricchezza_factor_w1 ricchezza_factor_w2 fisico_w1;
	output out = finale.all_stdz student = rmn2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w4 = ricchezza_factor_w1 ricchezza_factor_w2 fisico_w1 fisico_w2;
	output out = finale.all_stdz student = rn3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w4 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 fisico_w1 fisico_w2;
	output out = finale.all_stdz student = rmn3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 fisico_w1 fisico_w2 fisico_w4;
	output out = finale.all_stdz student = rn4;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 ricchezza_factor_w5 fisico_w1 fisico_w2 fisico_w4;
	output out = finale.all_stdz student = rmn4;
run;

* Denominator Weights for the Exposure and Mediator in the Outcome MSM ;
proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w2 = ricchezza_factor_w1 fisico_w1 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 country_w2 mstat_w2 age_w2 gender_w1 isced_w1;
	output out = finale.all_stdz student = rd2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w2 = ricchezza_factor_w1 ricchezza_factor_w2 fisico_w1 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 country_w2 mstat_w2 age_w2 bmi_w2 chronic_w2 gender_w1 isced_w1;
	output out = finale.all_stdz student = rmd2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w4 = ricchezza_factor_w1 ricchezza_factor_w2 fisico_w1 fisico_w2 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 country_w4 mstat_w4 age_w4 gender_w1 isced_w1;
	output out = finale.all_stdz student = rd3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w4 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 fisico_w1 fisico_w2 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 country_w4 mstat_w4 age_w4 bmi_w4 chronic_w4 gender_w1 isced_w1;
	output out = finale.all_stdz student = rmd3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 fisico_w1 fisico_w2 fisico_w4 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 benessere_w4 mstat_w4 country_w4 otrf1_w4 br015_w4 br016_w4 chronic_w4 bmi_w4 age_w4 country_w5 mstat_w5 age_w5 gender_w1 isced_w1;
	output out = finale.all_stdz student = rd4;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model fisico_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 ricchezza_factor_w5 fisico_w1 fisico_w2 fisico_w4 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 benessere_w4 mstat_w4 country_w4 otrf1_w4 br015_w4 br016_w4 chronic_w4 bmi_w4 age_w4 country_w5 mstat_w5 age_w5 bmi_w5 chronic_w5 gender_w1 isced_w1;
	output out = finale.all_stdz student = rmd4;
run;

* Form the Weights for the Exposure and Mediator in the Outcome MSM ;
data finale.all_stdz;
	set finale. all_stdz;
	wn2=2.506/(2.718**(-.5*rn2*rn2));
	wn3=2.506/(2.718**(-.5*rn3*rn3));
	wn4=2.506/(2.718**(-.5*rn4*rn4));
	wd2=2.506/(2.718**(-.5*rd2*rd2));
	wd3=2.506/(2.718**(-.5*rd3*rd3));
	wd4=2.506/(2.718**(-.5*rd4*rd4));
	wmn2=2.506/(2.718**(-.5*rmn2*rmn2));
	wmn3=2.506/(2.718**(-.5*rmn3*rmn3));
	wmn4=2.506/(2.718**(-.5*rmn4*rmn4));
	wmd2=2.506/(2.718**(-.5*rmd2*rmd2));
	wmd3=2.506/(2.718**(-.5*rmd3*rmd3));
	wmd4=2.506/(2.718**(-.5*rmd4*rmd4));
run;

data finale.all_stdz;
	set finale. all_stdz;
	ses = ricchezza_factor_w2 + ricchezza_factor_w4 + ricchezza_factor_w5;
	fisico = fisico_w2 + fisico_w4 + fisico_w5;
	weight = wd2 * wd3 * wd4 * wmd2 * wmd3 * wmd4 / (wn2 * wn3 * wn4 * wmn2 * wmn3 * wmn4);
	benessere = benessere_w6;
run;

* Outcome Marginal Structural Model for SES and Health on Well-being ;
ods output ParameterEstimates = parms1;
proc reg data = finale.all_stdz;
	model benessere = ses fisico / clb;
	weight weight;
run;
ods output close;

* MSM M ;
* Weights for Exposure History for the Mediator Marginal Structural Model ;
proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w2 = ricchezza_factor_w1;
	output out = finale.all_stdz student = rin2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w4 = ricchezza_factor_w1 ricchezza_factor_w2;
	output out = finale.all_stdz student = rin3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4;
	output out = finale.all_stdz student = rin4;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w2 = ricchezza_factor_w1 fisico_w1 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 country_w2 mstat_w2 age_w2 gender_w1 isced_w1;
	output out = finale.all_stdz student = rid2;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w4 = ricchezza_factor_w1 ricchezza_factor_w2 fisico_w1 fisico_w2 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 country_w4 mstat_w4 age_w4 gender_w1 isced_w1;
	output out = finale.all_stdz student = rid3;
run;

proc reg data = finale.all_stdz plots(maxpoints = none);
	model ricchezza_factor_w5 = ricchezza_factor_w1 ricchezza_factor_w2 ricchezza_factor_w4 fisico_w1 fisico_w2 fisico_w4 benessere_w1 mstat_w1 country_w1 otrf1_w1 br015_w1 br016_w1 chronic_w1 bmi_w1 age_w1 benessere_w2 mstat_w2 country_w2 otrf1_w2 br015_w2 br016_w2 chronic_w2 bmi_w2 age_w2 benessere_w4 mstat_w4 country_w4 otrf1_w4 br015_w4 br016_w4 chronic_w4 bmi_w4 age_w4 country_w5 mstat_w5 age_w5 gender_w1 isced_w1;
	output out = finale.all_stdz student = rid4;
run;

data finale.all_stdz;
	set finale.all_stdz;
	win2=2.506/(2.718**(-.5*rin2*rin2));
	win3=2.506/(2.718**(-.5*rin3*rin3));
	win4=2.506/(2.718**(-.5*rin4*rin4));
	wid2=2.506/(2.718**(-.5*rid2*rid2));
	wid3=2.506/(2.718**(-.5*rid3*rid3));
	wid4=2.506/(2.718**(-.5*rid4*rid4));
run;

data finale.all_stdz;
	set finale.all_stdz;
	ses2 = (ricchezza_factor_w2 + ricchezza_factor_w4 + ricchezza_factor_w5) / 3;
	weight2 = wid2 * wid3 * wid4 / (win2 * win3 * win4);
	fisico2 = fisico_w6;
run;

* Mediator Marginal Structural Model for Health on SES ;
ods output ParameterEstimates = parms2;
proc reg data = finale.all_stdz;
	model fisico2 = ses2 / clb;
	weight weight2;
run;
ods output close;

* Estimate Direct and Indirect Effects ;
data finale.temptt2;
	set parms1 parms2;
	bya = lag3(estimate);
	bym = lag3(estimate);
	direct = 3*lag3(estimate);
	indirect = 3*estimate*lag2(estimate);
	total = direct + indirect;
	keep direct indirect total;
	if _n_=5;
run; 

[FreelanceReinh · Accepted solution]

Hello @saramanzella,

My favorite paper on bootstrapping with SAS is David L. Cassell's BootstrapMania!: Re-Sampling the SAS® Way.

The key point is to use BY-processing rather than a macro loop. So, you would create a large input dataset containing all the bootstrap samples (or a view if a dataset would be too large, see p. 6 of the paper), insert a BY statement into each PROC REG step (and suppress the plots) and adapt the final DATA step correspondingly (as far as I see).

I would start with a few bootstrap samples to test the code and adjust settings (e.g. ods noresults) to optimize efficiency. Once everything is ready, increase the number of bootstrap samples to the final value.

Good luck with your thesis.