* 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:\sismid2020";

/* some preprocessing of the data */
data nhefs;
	set causinf.nhefs_book;
	if .<school<9 then education='1. 8th grade or less';
	else if 8<school<12 then education='2. HS dropout';
	else if school=12 then education='3. HS';
	else if 12<school<16 then education='4. College dropout';
	else if school ge 16 then education='5. College or more';
	else education= 'Unknown';
	university=	(school ge 16);
	older= (age>50);
	match=1;
	cens= (wt82 eq .);

	if education ne 'Unknown';
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;
	ods exclude ClassLevelInfo  Association FitStatistics GlobalTests;
	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 univariate data=nhefs_w;
	id seqn;
	var w;
run;

proc genmod data= nhefs_w; 
	class seqn;
	weight w;
	model wt82_71= qsmk;
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?

Sensitivity analysis: Four individuals have weights > 10; how does removing these individuals change the results?

*/