data  nhefs;

filename nhfile url
'http://publicifsv.sund.ku.dk/~pka/epidata/nhefs.txt';
infile nhfile firstobs=2;
input id wt82_71 qsmk sex race age education smokeintensity smokeyrs exercise active wt71 ht;
run;

/*
- wt82_71 (kgs)
- QUIT SMOKING BETWEEN 1ST QUESTIONNAIRE AND 1982, 1:YES, 0:NO
- 0: MALE 1: FEMALE
- 0: WHITE 1: BLACK OR OTHER IN 1971
- age (years)
- LEVEL OF EDUCATION IN 1971, 1: 8TH GRADE OR LESS,
  2: HS DROPOUT; 3: HS, 4: COLLEGE DROPOUT, 5: COLLEGE OR MORE
- NUMBER OF CIGARETTES SMOKED PER DAY IN 1971
- years smoked
- IN RECREATION, HOW MUCH EXERCISE? IN 1971,
  0:much exercise,1:moderate exercise,2:little or no exercise
- IN YOUR USUAL DAY, HOW ACTIVE ARE YOU? IN 1971,
  0:very active, 1:moderately active, 2:inactive, 3:missing
- wt71 (kgs)
- ht height (cms)
*/

/* 1: using the g-formula */

/* First we triple the data set */

data triple; set nhefs;
interv=-1; output;
interv=0; wt82_71=.; qsmk=0; output;
interv=1; wt82_71=.; qsmk=1; output;
run;

/* Next, we fit the Q-model - here a simple additive model */

proc genmod data=triple;
class education exercise active;
model wt82_71=qsmk sex race age education smokeintensity smokeyrs exercise
              active wt71 ht;
output out=preddata predicted=yhat;
run;

/* Finally, the mean outcome is predicted for qsmk=0 and qsmk=1 */

proc means data=preddata;
class interv;
var yhat;
run;

/* 2:  using IPTW */

/* First, we fit the propensity score model (simple additive) */

proc genmod data=nhefs descending;
class education exercise active;
model qsmk = sex race age education smokeintensity smokeyrs exercise
              active wt71 ht/dist=bin;
output out=iptw predicted=propscore;
run;

/* Next, we compute weights */

data iptw; set iptw;
if qsmk=1 then w=1/propscore;
if qsmk=0 then w=1/(1-propscore);
run;

/* and fit the marginal structural model using the weights */

proc genmod data=iptw;
class id;
model wt82_71 = qsmk;
repeated subject=id;
weight w;
run;