/**********************************************************************************************
 Program Name: PowerGxE.sas
 Program Date: 19 March 2007
 Programmer: Jean Yu-Chun Yen

%powerge calculates power in case-control studies examining genetic factors. 
It allows for power calculations for the marginal effects of gene,
as well as the joint test for genetic marginal effect and gene-environment interaction. 
Environmental effect is modelled as dichotomous. 
Genetic model is restricted to dominant inheritance model.

True penetrance model: 
logit[Pr(D=1|G,E)]=logit[base]+bg*G+be*E+bge*GE
G=1 for carriers of the risk allele and G=0 for non-carriers
E=1 for exposed subjects and E=0 for unexposed subjects

input   n:        Number of cases                  
        m:        Control:Case ratio               
        pg:       Risk allele frequency            
        pe:       Exposure frequency               
        base:     Baseline disease prevalence      
        bg:       "bg" in the true penetrance model  
        be:       "be" in the true penetrance model  
        bge:      "bge" in the true penetrance model 
        alpha:    Alpha level for the power calculation
        
output  POWER data set
        Pow_GGE:  Power for marginal genetic test
        Pow_G:    Power for joint test
        Pow_GE:   Power for gene-environment interaciton test
        pg:       Risk allele frequency
        pe:       Exposure frequency
        Rg:       exp(bg) in the true penetrance model
        Re:       exp(be) in the true penetrance model
        Rge:      exp(bge) in the true penetrance model 
        alpha:    Alpha level for the power calculation
**********************************************************************************************/

/*******************************************************************/
/*                                DEFINE MACROFUNCTION POWERGE     */
%macro powerge(n=,             /* Number of cases                  */
               m=1,            /* Control:Case ratio               */
               pg=,            /* Risk allele frequency            */
               pe=,            /* Exposure frequency               */
               base=,          /* baseline disease prevalence      */
               bg=,            /* bg in the true penetrance model  */
               be=,            /* be in the true penetrance model  */
               bge=,           /* bge in the true penetrance model */
               alpha=.05);     /* alpha level                      */
/*******************************************************************/

data test;
  p=&pg; q=1-p;
  alpha=log(&base/(1-&base)); bg=&bg; be=&be; bge=&bge;
  array Z{3}   _temporary_ (0 1 1);  *Genotype coding for dominant inheritance model;
  array PrG{3} _temporary_;
  PrG{1}=q**2; PrG{2}=p*q*2; PrG{3}=p**2;
  array PrE{2} _temporary_;
  PrE{1}=1-&pe; PrE{2}=&pe;
  * prevalence;
  totprob=0;
  do e=0 to 1;
  do g=1 to 3;
	 totprob=totprob+PrE{e+1}*PrG{g}*exp(alpha+bg*Z{g}+be*e+bge*Z{g}*e)/(1+exp(alpha+bg*Z{g}+be*e+bge*Z{g}*e));
  end;
  end;
  * cases;
  y=1;
  do e=0 to 1;
  do g=1 to 3;
     PrDG=PrG{g}*PrE{e+1}*exp(alpha+bg*Z{g}+be*e+bge*Z{g}*e)/(1+exp(alpha+bg*Z{g}+be*e+bge*Z{g}*e));
     PrDG=PrDG/totprob;
	 zg=Z{g};
	 output;
  end;
  end;
  * controls;
  y=0;
  do e=0 to 1;
  do g=1 to 3;
     PrDG=PrG{g}*PrE{e+1}/(1+exp(alpha+bg*Z{g}+be*e+bge*Z{g}*e));
     PrDG=&m*PrDG/(1-totprob);
	 zg=Z{g};
	 output;
  end;
  end;
  keep y zg e PrDG;
run;

proc logistic noprint data=test outest=L0;
   model y(descending)=;
   weight PrDG;
proc logistic noprint data=test outest=L1;
   model y(descending)=e zg e*zg;
   weight PrDG;
proc logistic noprint data=test outest=LE;
   model y(descending)=e;
   weight PrDG;
proc logistic noprint data=test outest=LGE;
   model y(descending)=e zg;
   weight PrDG;
proc logistic noprint data=test outest=LG;
   model y(descending)=zg;
   weight PrDG;


data POWER;
   merge L0(rename=(_LNLIKE_=like0)) L1(rename=(_LNLIKE_=like1)) LE(rename=(_LNLIKE_=likee))
         LG(rename=(_LNLIKE_=likeg)) LGE(rename=(_LNLIKE_=likege));
   
   *G-GE;
   NCP_GGE=&n*2*(like1-likee);
   Pow_GGE=1-probchi(cinv(1-&alpha,2),2,NCP_GGE);
   
   *G;
    NCP_G=&n*2*(likeg-like0);
    Pow_G=1-probchi(cinv(1-&alpha,1),1,NCP_G);
    
   *GE;
    NCP_GE=&n*2*(like1-likege);
    Pow_GE=1-probchi(cinv(1-&alpha,1),1,NCP_GE);
   
   pg=&pg; Re=exp(&be); Rg=exp(&bg); Rge=exp(&bge); pe=&pe; alpha=α n=&n; m=&m;
   keep Pow_GGE Pow_G Pow_GE pg pe Rg Re Rge alpha n m;
run;

proc datasets; delete test L0 L1 LE LGE LG;

ods select all;

proc print data=POWER noobs;
  title "Power for case-control studies";
run;
%mend;

/* Example */
/* 
%powerge(n=500,m=1,pg=.1,pe=.25,base=.0001,bg=log(1.0),be=log(1.5),bge=log(1.0),alpha=.05);
%powerge(n=500,m=1,pg=.1,pe=.25,base=.0001,bg=log(1.0),be=log(1.5),bge=log(2.0),alpha=.05);
*/