* Macro MOFET fits case-parent model described in Kraft et al. 2003 * Peter Kraft * 7/14/2003: First Draft * 7/14/2004: Fixed maternal-effect weights (multiple cases), denominator and numerator * 7/29/2004: Fixed interaction-effect weights (multiple cases), denominator and numerator; %macro mofet(data=in, d=case, g=geno, momg=momg, dadg=dadg, mt=mt, nfam=nfam, maxnfam=7, parms=offmain mfint = 0, model=offmain*Off + mfint*MomD*Off, nlmopt= ); data UseMe; set &data; d=1; run; * Initialize dummy outcome; proc nlmixed data=UseMe &nlmopt; array Gee{&maxnfam} &g.1-&g.&maxnfam; array Geed{&maxnfam} &g.d1-&g.d&maxnfam; array Dee{&maxnfam} &d.1-&d.&maxnfam; * array Eee{7} Exposure1-Exposure7; array PrG_MT0{3} _temporary_; PrG_MT0{1} = 1; PrG_MT0{2} = 0; PrG_MT0{3} = 0; array PrG_MT2{3} _temporary_; PrG_MT2{1} =.5; PrG_MT2{2} =.5; PrG_MT2{3} = 0; array PrG_MT4{3} _temporary_; PrG_MT4{1} = 0; PrG_MT4{2} = 1; PrG_MT4{3} = 0; array PrG_MT5{3} _temporary_; PrG_MT5{1} =.25; PrG_MT5{2} =.5; PrG_MT5{3} =.25; array PrG_MT6{3} _temporary_; PrG_MT6{1} = 0; PrG_MT6{2} =.5; PrG_MT6{3} =.5; array PrG_MT8{3} _temporary_; PrG_MT8{1} = 0; PrG_MT8{2} = 0; PrG_MT8{3} = 1; exptheta0=1-exp(theta2)-exp(theta4)-exp(theta5)-exp(theta6)-exp(theta8); PrMT0=exptheta0; PrMT2=exp(theta2); PrMT4=exp(theta4); PrMT5=exp(theta5); PrMT6=exp(theta6); PrMT8=exp(theta8); ncases=sum(of &d.1-&d.&maxnfam); /* BEGIN define Likelihood denom, Pr(D) */ denom=0; do m=0 to 2; do f=0 to 2; * loop over all possible parental genotypes; /* initialize maternal genotype scores, mating type */ MomD=(m ge 1)*ncases; MomR=(m eq 2)*ncases; Mom=m*ncases; Mom1=(m eq 1)*ncases; Mom2=(m eq 2)*ncases; mtd=2*(m+f) + ((m eq 1) and (f eq 1)); do i=0 to ((3**&nfam)-1); * loop over all possible offspring genotypes; /* initialize offspring genotypes, genotype scores, calculate Pr(Go|Gp) */ OffD=0; OffR=0; Off=0; Off1=0; Off2=0; PrGee_MT0=1; PrGee_MT2=1; PrGee_MT4=1; PrGee_MT5=1; PrGee_MT6=1; PrGee_MT8=1; temp=i; do j=1 to &nfam; Geed{j}=mod(temp,3); OffD=OffD+(Geed{j} ge 1)*Dee{j}; OffR=OffR+(Geed{j} eq 2)*Dee{j}; Off =Off +Geed{j}*Dee{j}; Off1=Off1+(Geed{j} eq 1)*Dee{j}; Off2=Off2+(Geed{j} eq 2)*Dee{j}; IntD=OffD*(m ge 1); IntR=OffD*(m eq 2); Int =Off * m; Int11=Off1*(m eq 1); Int12=Off1*(m eq 2); Int21=Off2*(m eq 1); Int22=Off2*(m eq 2); PrGee_MT0=PrGee_MT0*PrG_MT0{Geed{j}+1}; PrGee_MT2=PrGee_MT2*PrG_MT2{Geed{j}+1}; PrGee_MT4=PrGee_MT4*PrG_MT4{Geed{j}+1}; PrGee_MT5=PrGee_MT5*PrG_MT5{Geed{j}+1}; PrGee_MT6=PrGee_MT6*PrG_MT6{Geed{j}+1}; PrGee_MT8=PrGee_MT8*PrG_MT8{Geed{j}+1}; temp=floor(temp/3); end; /* calculate prob */ /* terms in denominator summations have form Pr(MT) Pr(Gp|MT) Pr(Go|Gp) Pr(Dcase|Go,Gp) */ if (mtd eq 0) then Prob = PrMT0 *PrGee_MT0; if (mtd eq 2) then Prob = PrMT2*.5*PrGee_MT2*exp(&model); if (mtd eq 4) then Prob = PrMT4*.5*PrGee_MT4*exp(&model); if (mtd eq 5) then Prob = PrMT5 *PrGee_MT5*exp(&model); if (mtd eq 6) then Prob = PrMT6*.5*PrGee_MT6*exp(&model); if (mtd eq 8) then Prob = PrMT8 *PrGee_MT8*exp(&model); denom=denom+Prob; end; end; end; /* END define Likelihood denom, Pr(D) */ /* BEGIN define Likelihood num, Pr(D|Go,Gm) Pr(Go|Gp) Pr(Gp) */ MomD=(&momg ge 1)*ncases; MomR=(&momg eq 2)*ncases; Mom=&momg.*ncases; Mom1=(&momg eq 1)*ncases; Mom2=(&momg eq 2)*ncases; OffD=0; OffR=0; Off=0; Off1=0; Off2=0; PrGee_MT0=1; PrGee_MT2=1; PrGee_MT4=1; PrGee_MT5=1; PrGee_MT6=1; PrGee_MT8=1; do i=1 to nfam; OffD=OffD+(Gee{i} ge 1)*Dee{i}; OffR=OffR+(Gee{i} eq 2)*Dee{i}; Off =Off +Gee{i}*Dee{i}; Off1=Off1+(Gee{i} eq 1)*Dee{i}; Off2=Off2+(Gee{i} eq 2)*Dee{i}; IntD=OffD*(&momg ge 1); IntR=OffD*(&momg eq 2); Int =Off * &momg; Int11=Off1*(&momg eq 1); Int12=Off1*(&momg eq 2); Int21=Off2*(&momg eq 1); Int22=Off2*(&momg eq 2); PrGee_MT0=PrGee_MT0*PrG_MT0{Gee{i}+1}; PrGee_MT2=PrGee_MT2*PrG_MT2{Gee{i}+1}; PrGee_MT4=PrGee_MT4*PrG_MT4{Gee{i}+1}; PrGee_MT5=PrGee_MT5*PrG_MT5{Gee{i}+1}; PrGee_MT6=PrGee_MT6*PrG_MT6{Gee{i}+1}; PrGee_MT8=PrGee_MT8*PrG_MT8{Gee{i}+1}; end; if (&mt eq 0) then num = PrMT0 *PrGee_MT0; if (&mt eq 2) then num = PrMT2*.5*PrGee_MT2*exp(&model); if (&mt eq 4) then num = PrMT4*.5*PrGee_MT4*exp(&model); if (&mt eq 5) then num = PrMT5 *PrGee_MT5*exp(&model); if (&mt eq 6) then num = PrMT6*.5*PrGee_MT6*exp(&model); if (&mt eq 8) then num = PrMT8 *PrGee_MT8*exp(&model); /* BEGIN MISSSING FATHER families */ if (&mt eq 10) then do; if (&momg eq 0) then do; num= PrMT0 *PrGee_MT0; num=num+PrMT2*.5*PrGee_MT2*exp(&model); num=num+PrMT4*.5*PrGee_MT4*exp(&model); end; if (&momg eq 1) then do; num= PrMT2*.5*PrGee_MT2*exp(&model); num=num+PrMT5 *PrGee_MT5*exp(&model); num=num+PrMT6*.5*PrGee_MT6*exp(&model); end; if (&momg eq 2) then do; num= PrMT4*.5*PrGee_MT4*exp(&model); num=num+PrMT6*.5*PrGee_MT6*exp(&model); num=num+PrMT8 *PrGee_MT8*exp(&model); end; end; /* END MISSSING FATHER families */ /* BEGIN MISSSING MOTHER families */ if (mt eq 12) then do; if (&dadg eq 0) then do; num= PrMT0 *PrGee_MT0; MomD=1*ncases; MomR=0*ncases; Mom=1*ncases; Mom1=1*ncases; Mom2=0*ncases; num=num+PrMT2*.5*PrGee_MT2*exp(&model); MomD=1*ncases; MomR=1*ncases; Mom=2*ncases; Mom1=0*ncases; Mom2=1*ncases; num=num+PrMT4*.5*PrGee_MT4*exp(&model); end; if (&dadg eq 1) then do; MomD=0*ncases; MomR=0*ncases; Mom=0*ncases; Mom1=0*ncases; Mom2=0*ncases; num= PrMT2*.5*PrGee_MT2*exp(&model); MomD=1*ncases; MomR=0*ncases; Mom=1*ncases; Mom1=1*ncases; Mom2=0*ncases; num=num+PrMT5 *PrGee_MT5*exp(&model); MomD=1*ncases; MomR=1*ncases; Mom=2*ncases; Mom1=0*ncases; Mom2=1*ncases; num=num+PrMT6*.5*PrGee_MT6*exp(&model); end; if (&dadg eq 2) then do; MomD=0*ncases; MomR=0*ncases; Mom=0*ncases; Mom1=0*ncases; Mom2=0*ncases; num= PrMT4*.5*PrGee_MT4*exp(&model); MomD=1*ncases; MomR=0*ncases; Mom=1*ncases; Mom1=1*ncases; Mom2=0*ncases; num=num+PrMT6*.5*PrGee_MT6*exp(&model); MomD=1*ncases; MomR=1*ncases; Mom=2*ncases; Mom1=0*ncases; Mom2=1*ncases; num=num+PrMT8 *PrGee_MT8*exp(&model); end; end; /* END MISSSING MOTHER families */ /* END define Likelihood num, Pr(D|Go,Gm) Pr(Go|Gp) Pr(Gp) */ ll = d*(log(num)-log(denom)); model d ~ general(ll); parms theta2 theta4 theta5 theta6 theta8 = -1.7918 ; ods output parameterestimates=parm; run; %mend;