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

/***************************************************************
PROGRAM 14.1
Preprocessing, ranks of extreme observations, IP weights for censoring
Data from NHEFS
***************************************************************/;
libname causinf "V:\sismid2023";

/* some preprocessing of the data */
data nhefs;
    set causinf.nhefs;
    cens= (wt82 eq .);
run;

/* estimation of denominator of ip weights for C */
proc logistic data=nhefs;
    ods exclude ClassLevelInfo  Association FitStatistics GlobalTests;
    class exercise active education;
    model cens = qsmk sex race age age*age education
                smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs
                exercise active wt71 wt71*wt71;
    output out=est_prob_d_c (keep= seqn pd_cens) p=pd_cens;
run;
proc sort; by seqn; run;

data nhefs_w;
    merge nhefs est_prob_d_c;
    by seqn;
    if cens=0; * observations with cens=1 only contribute to censoring models;
    w_c= 1 / pd_cens;
run;

proc univariate data=nhefs_w;
    var w_c;
    id seqn;
run;


/*****************************************************************
PROGRAM 14.2
G-estimation of a 1-parameter structural nested mean model
Brute force search
Data from NHEFS
******************************************************************/;

/*********************************************************/
/* G-estimation: Checking multiple possible values of psi*/
/*********************************************************/

%macro gest_snm;

    %do counter=200 %to 500 %by 10; /* integers only in do loops */
        %let psi= %eval(&counter)/100;
        data nhefs_gest;
            set nhefs_w;
            Hpsi= wt82_71 - &psi*qsmk;
        run;

        title "Testing PSI= &psi";

        proc genmod data= nhefs_gest descending; 
            ods exclude ClassLevels ParmInfo ModelInfo;  
            class seqn exercise active education;
            weight w_c;
            model qsmk = sex race age age*age education
                        smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs
                        exercise active wt71 wt71*wt71
                        Hpsi/ dist= bin link= logit;
            repeated subject=seqn / type=ind;
        run;
        quit;

    %end;

%mend;

%gest_snm;
title;

/*********************************************************/
/* G-estimation: Checking one possible value of psi      */
/*********************************************************/
data nhefs_gest;
    set nhefs_w;
    psi= 3.446;
    Hpsi= wt82_71-psi*qsmk;
run;

run;

proc genmod data= nhefs_gest descending;
    ods exclude ClassLevels ParmInfo ModelInfo;  
    class seqn exercise active education;
    weight w_c;
    model qsmk = sex race age age*age education
                smokeintensity smokeintensity*smokeintensity smokeyrs smokeyrs*smokeyrs
                exercise active wt71 wt71*wt71
                Hpsi/ dist= bin link= logit;
    repeated subject=seqn / type=ind;
run;
quit;