* Adapted from programs found here https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/ ;

/***************************************************************
PROGRAM 12.1
Descriptive statistics from NHEFS data
***************************************************************/;
libname causinf "V:\sismid2023";

/* some preprocessing of the data */
data nhefs;
	set causinf.nhefs;
run;

data nhefs_nmv;
	set nhefs;
	if wt82 ne .; * provisionally ignore subjects with missing values for weight in 1982;
run;

proc genmod data= nhefs_nmv; 
	model wt82_71= qsmk;
	estimate 'Smoking cessation' intercept 1 qsmk 1;
	estimate 'No smoking cessation' intercept 1 qsmk 0;
run;



/***************************************************************
PROGRAM 12.2
Estimating IP weights
Data from NHEFS
***************************************************************/;

/* estimation of ip weights via a logistic model */
proc logistic data= nhefs_nmv descending;
	class exercise active education;
	model qsmk = sex race age age*age education
				smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs
				exercise active wt71 wt71*wt71;
	output out=est_prob p=p_qsmk;
run;

data nhefs_w;
	set est_prob;
	if qsmk=1 then w= 1/p_qsmk;
	else if qsmk=0 then w= 1/(1-p_qsmk);
	
run;


proc genmod data= nhefs_w; 
	class seqn;
	weight w;
	model wt82_71= qsmk;
	* repeated subject=seqn / type=ind;
run;

/*

What is IPW estimate of effect of qsmk?

What is IPW estimate of mean outcome if all individuals quit smoking? If no individuals quit smoking?

*/