*************************************************************************************************;
*        CLIENT: Moderna P301 
*      PROTOCOL: mRNA-1273-P301
*       PURPOSE: Create Inferential Statistical Analysis Macro for Primary and Secondary Efficacy Analyses
*   INPUT FILES: ADaM datasets 
*  OUTPUT FILES: stat macro
*   USAGE NOTES:  
*************************************************************************************************;
*  Copyright 2020 PPD
*  All Rights Reserved. 
*************************************************************************************************;


%macro VEHRIR(inlib=adb,table=,dataset=,subset=,treatment=,time=,case=);
data vehr_&table.; run;
data veir_&table.; run;

proc sql noprint;
    create table _temp1 as
    select a.*, b.BMIBL
    from &inlib..&dataset. as a left join adb.adsl as b
    on a.USUBJID=b.USUBJID
    ;
quit;

data temp1;
set _temp1;
if &subset.;
run;


/* run only if obs gt 0*/
     proc sql noprint;
      select count(*) 
        into :obs_count
          from temp1;
    quit;

    %if &obs_count = 0 %then
      %do;
	  %put No observation;
	%end;

    %if &obs_count gt 0 %then
      %do;

%let table2=%index(%upcase(&subset.),STRATAVN)+%index(%upcase(&subset.),AGEGR1)*%index(%upcase(&subset.),RISKGR1);
%put &table2;
%if &table2=0 %then 
  %do;
proc phreg data=temp1;
class &treatment. (ref='Placebo');
model &time.*&case.(1)= &treatment./alpha=0.05 ties=Efron rl;
estimate &treatment. 1/testvalue=-0.3566749 lower cl;
strata STRATAR;
ods output ParameterEstimates=HR Estimates=Pval(keep=Probz);
run;
%end;
%else %if &table2.>0 %then 
   %do;
proc phreg data=temp1;
class &treatment. (ref='Placebo');
model &time.*&case.(1)= &treatment./alpha=0.05 ties=Efron rl;
estimate &treatment. 1/testvalue=-0.3566749 lower cl;
ods output ParameterEstimates=HR Estimates=Pval(keep=Probz);
run;
%end;

/*data infout.VEHR_&table.;*/
data VEHR_&table.;
merge HR PVAL;
VE=1-HazardRatio;
VELOW= 1-HRUpperCL;
VEUPP= 1-HRLowerCL;
VEPVAL=Probz;
run;

/*incidence rate*/
/*number of cases*/
proc freq data=temp1;
table &treatment.*&case.;
ods output CrossTabFreqs=freqs;
run;
data freqs2;
set freqs;
if not missing(&treatment.) and not missing(rowpercent);
if rowpercent=100 and &case.=1 then 
  do;
&case._0=0; Frequency_0=0;rowpercent_0=0;
drop &case. Frequency RowPercent;
rename &case._0=CNSR Frequency_0=Frequency RowPercent_0=RowPercent;
end;
else 
  do;
&case._0=&case.; Frequency_0=Frequency;rowpercent_0=rowpercent;
drop &case. Frequency RowPercent;
rename &case._0=CNSR Frequency_0=Frequency RowPercent_0=RowPercent;
end;
run;
data freqs2;
set freqs2;
keep &treatment. &case. Frequency RowPercent;
run;
proc sort data=freqs2 nodupkey;
by &treatment. &case. Frequency RowPercent;
run;
/* person year*/
proc means data=temp1 sum;
var &time.;
class &treatment.;
ods output summary=sumpm;
run;
data freqs3;
set freqs2;
if &case.=0;
run;

/* Prepare GENMOD for IR dataset*/
proc sort data=freqs3;
by &treatment.;
run;
proc sort data=sumpm;
by &treatment.;
run;

data IRdata;
merge freqs3 sumpm;
by &treatment.;
logpm=log(&time._SUM/365.25);
run;

/*Cruide IR*/
data IROUT;
set IRdata;
pm=&time._sum/365.25;
if frequency>0 then caseLow=cinv(.025, (frequency*2))/2;
caseupp=cinv(.975, ((frequency+1)*2))/2;
IR=Frequency/pm*1000;
IRLOW=caselow/pm*1000;
IRUPP=caseupp/pm*1000;
run;

/*Calcualte the exact CI for IRR*/
data IRDATA_P;
set IRDATA;
if &treatment.='Placebo';
casep=frequency;
sumpmP=&time._sum/365.25;
keep caseP sumpmP;
run;

data IRDATA_M;
set IRDATA;
if &treatment.='mRNA-1273';
caseM=frequency;
sumpmM=&time._sum/365.25;
keep caseM sumpmM;
run;
data IRdataCI;
merge IRDATA_P IRDATA_M;
if caseP>0 then IRR=(caseM/sumpmM)/(caseP/sumpmP);
if caseM>0 then f_low=FINV(0.975,(caseP+1)*2,caseM*2);
if caseP>0 then f_upp=FINV(0.975,(caseM+1)*2,caseP*2);
if caseP>0 then IRR_E=(caseM/sumpmM)/(caseP/sumpmP);
if caseP>0 then IRRLOW_E=(sumpmP/sumpmM)*(caseM/(caseP+1))*(1/f_low);
if caseP>0 then IRRUPP_E=(sumpmP/sumpmM)*((caseM+1)/caseP)*f_upp;
&treatment.='mRNA-1273';
keep &treatment. IRR_E IRRLOW_E IRRUPP_E;
run;

proc sort data=IROUT;
by &treatment.;
run;
proc sort data=IRDATACI;
by &treatment.;
run;
/*data infout.VEIR_&table.;*/
data VEIR_&table.;
merge IROUT IRDATACI;
by &treatment.;
VEIR=1-IRR_E;
VEIRLOW=1-IRRUPP_E;
VEIRUPP=1-IRRLOW_E;
run;
%end;
%mend;



%macro VEHRIR2(inlib=adb,table=,dataset=,subset=,treatment=,time=,case=,strata=);
data vehr_&table.; run;
data veir_&table.; run;
 
data temp1;
set &inlib..&dataset.;
if &subset.;
run;


/* run only if obs gt 0*/
     proc sql noprint;
      select count(*) 
        into :obs_count
          from temp1;
    quit;

    %if &obs_count = 0 %then
      %do;
	  %put No observation;
	%end;

    %if &obs_count gt 0 %then
      %do;

%let table2=%index(%upcase(&subset.),STRATAVN)+%index(%upcase(&subset.),AGEGR1)*%index(%upcase(&subset.),RISKGR1);
%put &table2;
%if &table2=0 %then 
  %do;
proc phreg data=temp1;
class &treatment. (ref='Placebo');
model &time.*&case.(1)= &treatment./alpha=0.05 ties=Efron rl;
estimate &treatment. 1/testvalue=-0.3566749 lower cl;
strata &strata.;
ods output ParameterEstimates=HR Estimates=Pval(keep=Probz);
run;
%end;
%else %if &table2.>0 %then 
   %do;
proc phreg data=temp1;
class &treatment. (ref='Placebo');
model &time.*&case.(1)= &treatment./alpha=0.05 ties=Efron rl;
estimate &treatment. 1/testvalue=-0.3566749 lower cl;
ods output ParameterEstimates=HR Estimates=Pval(keep=Probz);
run;
%end;

/*data infout.VEHR_&table.;*/
data VEHR_&table.;
merge HR PVAL;
VE=1-HazardRatio;
VELOW= 1-HRUpperCL;
VEUPP= 1-HRLowerCL;
VEPVAL=Probz;
run;

/*incidence rate*/
/*number of cases*/
proc freq data=temp1;
table &treatment.*&case.;
ods output CrossTabFreqs=freqs;
run;
data freqs2;
set freqs;
if not missing(&treatment.) and not missing(rowpercent);
if rowpercent=100 and &case.=1 then 
  do;
&case._0=0; Frequency_0=0;rowpercent_0=0;
drop &case. Frequency RowPercent;
rename &case._0=CNSR Frequency_0=Frequency RowPercent_0=RowPercent;
end;
else 
  do;
&case._0=&case.; Frequency_0=Frequency;rowpercent_0=rowpercent;
drop &case. Frequency RowPercent;
rename &case._0=CNSR Frequency_0=Frequency RowPercent_0=RowPercent;
end;
run;
data freqs2;
set freqs2;
keep &treatment. &case. Frequency RowPercent;
run;
proc sort data=freqs2 nodupkey;
by &treatment. &case. Frequency RowPercent;
run;
/* person year*/
proc means data=temp1 sum;
var &time.;
class &treatment.;
ods output summary=sumpm;
run;
data freqs3;
set freqs2;
if &case.=0;
run;

/* Prepare GENMOD for IR dataset*/
proc sort data=freqs3;
by &treatment.;
run;
proc sort data=sumpm;
by &treatment.;
run;

data IRdata;
merge freqs3 sumpm;
by &treatment.;
logpm=log(&time._SUM/365.25);
run;

/*Cruide IR*/
data IROUT;
set IRdata;
pm=&time._sum/365.25;
if frequency>0 then caseLow=cinv(.025, (frequency*2))/2;
caseupp=cinv(.975, ((frequency+1)*2))/2;
IR=Frequency/pm*1000;
IRLOW=caselow/pm*1000;
IRUPP=caseupp/pm*1000;
run;

/*Calcualte the exact CI for IRR*/
data IRDATA_P;
set IRDATA;
if &treatment.='Placebo';
casep=frequency;
sumpmP=&time._sum/365.25;
keep caseP sumpmP;
run;

data IRDATA_M;
set IRDATA;
if &treatment.='mRNA-1273';
caseM=frequency;
sumpmM=&time._sum/365.25;
keep caseM sumpmM;
run;
data IRdataCI;
merge IRDATA_P IRDATA_M;
if caseP>0 then IRR=(caseM/sumpmM)/(caseP/sumpmP);
if caseM>0 then f_low=FINV(0.975,(caseP+1)*2,caseM*2);
if caseP>0 then f_upp=FINV(0.975,(caseM+1)*2,caseP*2);
if caseP>0 then IRR_E=(caseM/sumpmM)/(caseP/sumpmP);
if caseP>0 then IRRLOW_E=(sumpmP/sumpmM)*(caseM/(caseP+1))*(1/f_low);
if caseP>0 then IRRUPP_E=(sumpmP/sumpmM)*((caseM+1)/caseP)*f_upp;
&treatment.='mRNA-1273';
keep &treatment. IRR_E IRRLOW_E IRRUPP_E;
run;

proc sort data=IROUT;
by &treatment.;
run;
proc sort data=IRDATACI;
by &treatment.;
run;
/*data infout.VEIR_&table.;*/
data VEIR_&table.;
merge IROUT IRDATACI;
by &treatment.;
VEIR=1-IRR_E;
VEIRLOW=1-IRRUPP_E;
VEIRUPP=1-IRRLOW_E;
run;
%end;
%mend;