SAS Macro for Testing Maternally-mediated Genetic Effects
/* READ COMMENTS:
This is a 2-degree-of-freedom test for maternally-mediated effects with missing fathers or mothers.  

Should be used with SAS Version 8 or higher. Runs in batch.  

This code works for the following SAS data.
Each record must consist of M,F,C, and X variables, where 
M= # of copies of the variant allele that the mother carries,
F= # of copies of the variant allele that the father carries,
C= # of copies of the variant allele that the child carries, and
X= the quantitative trait.
FLAG_F= flag that equals 0 if father's genotype is available and 1 if father's genotype is missing.
FLAG_M= flag that equals 0 if mother's genotype is available and 1 if mother's genotype is missing.  
ID= identifier of trio.

ENTER THE NAME OF YOUR DATA SET IN THE FIRST %LET STATEMENT BELOW.

If you want to increase the number of iterations, increase the variable K in the macro calls 

%emalg(X MTYPE01 MTYPE11 MTYPE12 MTYPE00 MTYPE02 MTYPE22,X,MTYPE01 MTYPE12 MTYPE02 INDX1 INDX2,k=25);

There are two calls to the macro in the body of the program.  Search on the number 25 and replace all with the new number of iterations. 
It is advisable to confirm convergence by trying a larger number of iterations, to be sure results are unchanged. 

If you want to test parent-of-origin effects change the third string of variables in the macro call, such that in the first call you model:
MTYPE01 MTYPE12 MTYPE02 INDX1 INDX2 K1X

while in the second call leave out the associated predictor K1X and model:
MTYPE01 MTYPE12 MTYPE02 INDX1 INDX2

You will need to change the degrees of freedom to 1 in the last data set for the line reading, ProbChiSq=1-PROBCHI(LRT,2);.  

Beware: There is no warranty with this code. 
*/

TITLE1 "POLYTOMOUS logistic approach for missing fathers or mothers";

%LET YOUR_DATA=ENTERHERE;
ODS LISTING CLOSE;
%MACRO emalg(vars1,vars2,maternal,k=);

%DO i=0 %TO &k;

%IF &i=0 %THEN %DO;
    %LET COND=((flag_f=0) and (flag_m=0));
    %LET D1=;
    %LET D2=;
    %END;
%ELSE %DO;
    %LET D1=B0 B2 C000 C002 C010 C012 C110 C112 C020 C022 C120 C122 C120 C122 C220 C222;
    %LET D2=A00 A01 A11 A12 A02 A22 U00 U01 U11 U12 U02 U22 mom1 mom2 imprint 
                                                            intmom01 intmom02 intmom12;
    %LET COND=((flag_f=0) OR (flag_m=0));
       PROC SORT DATA=WEIGHTS&j OUT=PLOGDATA&i;
       BY ORDER;
    %END;

*M-STEP;


%IF (&i=0) %THEN %DO;
PROC SORT DATA=PLOGDATA&i(where=(((((M=0) AND (F=2)) OR ((M=2) AND (F=0)) OR 
                                   ((M=0) AND (F=1)) OR ((M=1) AND (F=0)) OR 
                                   ((M=2) AND (F=1)) OR ((M=1) AND (F=2)))))) OUT=MDATA&i;
BY EFFECT;
%END;

%IF (&i^=0) %THEN %DO;
PROC SORT DATA=WEIGHTS&j(where=(((((M=0) AND (F=2)) OR ((M=2) AND (F=0)) OR 
                                  ((M=0) AND (F=1)) OR ((M=1) AND (F=0)) OR 
                                  ((M=2) AND (F=1)) OR ((M=1) AND (F=2)))))) OUT=MDATA&i;
BY EFFECT;
%END;


ODS TRACE OFF;
ODS OUTPUT ParameterEstimates=Mparms1ests&i;

PROC LOGISTIC DATA=MDATA&i(where=&COND) descending;
MODEL EFFECT = &maternal / NOINT link=logit;
WEIGHT W;
RUN; 

data Mparms&i;
  set mparms1ests&i END=EOF;
  RETAIN mom1 0 mom2 0 stemom1 0 stemom2 0 imprint 0 intmom01 0 intmom02 0 intmom12 0;
  IF VARIABLE="INDX1" THEN mom1=ESTIMATE;
  IF VARIABLE="INDX1" THEN stemom1=STDERR;
  IF VARIABLE="INDX2" THEN mom2=ESTIMATE;
  IF VARIABLE="INDX2" THEN stemom2=STDERR;
  IF VARIABLE="K1X" THEN imprint=ESTIMATE;
  IF VARIABLE="MTYPE01" THEN intmom01=ESTIMATE;
  IF VARIABLE="MTYPE12" THEN intmom12=ESTIMATE;
  IF VARIABLE="MTYPE02" THEN intmom02=ESTIMATE;
  IF EOF THEN OUTPUT;
  DROP ESTIMATE VARIABLE DF STDERR WALDCHISQ PROBCHISQ;

run;

ODS TRACE OFF;
ODS OUTPUT FitStatistics=lldata&i  ParameterEstimates=Parm1Ests&i;

PROC LOGISTIC DATA=PLOGDATA&i(where=&COND) ORDER=DATA;
MODEL C = &vars1 / NOINT LINK=GLOGIT;
WEIGHT W;
RUN;

DATA PARMS&i;
  SET Parm1Ests&i END=EOF;
  RETAIN B0 0 B2 0 C010 C012 C110 C112 C120 C122 C000 C002 C020 C022 C220 C222;
  IF VARIABLE="X" AND RESPONSE="0" THEN B0=ESTIMATE;
  IF VARIABLE="X" AND RESPONSE="2" THEN B2=ESTIMATE;
  IF VARIABLE="MTYPE01" AND RESPONSE="0" THEN C010=ESTIMATE;
  IF VARIABLE="MTYPE01" AND RESPONSE="2" THEN C012=ESTIMATE;
  IF VARIABLE="MTYPE11" AND RESPONSE="0" THEN C110=ESTIMATE;
  IF VARIABLE="MTYPE11" AND RESPONSE="2" THEN C112=ESTIMATE;
  IF VARIABLE="MTYPE12" AND RESPONSE="0" THEN C120=ESTIMATE;
  IF VARIABLE="MTYPE12" AND RESPONSE="2" THEN C122=ESTIMATE;
  IF VARIABLE="MTYPE00" AND RESPONSE="0" THEN C000=ESTIMATE;
  IF VARIABLE="MTYPE00" AND RESPONSE="2" THEN C002=ESTIMATE;
  IF VARIABLE="MTYPE02" AND RESPONSE="0" THEN C020=ESTIMATE;
  IF VARIABLE="MTYPE02" AND RESPONSE="2" THEN C022=ESTIMATE;
  IF VARIABLE="MTYPE22" AND RESPONSE="0" THEN C220=ESTIMATE;
  IF VARIABLE="MTYPE22" AND RESPONSE="2" THEN C222=ESTIMATE;

  IF EOF THEN OUTPUT;
  DROP ESTIMATE VARIABLE RESPONSE DF STDERR WALDCHISQ PROBCHISQ;
RUN;
/*
ods listing;
proc print data=parms&i;
var b0 b2;
run;
ods listing close;
*/
DATA CONDPROB&i;
  IF _N_=1 THEN SET PARMS&i;
  ELSE SET PLOGDATA&i(DROP=&D1) END=EOF;
  IF _N_^=1 THEN DO;
     P0_00=EXP(B0*X+C000)/(EXP(B0*X+C000)+1+EXP(B2*X+C002));
     P0_01=EXP(B0*X+C010)/(EXP(B0*X+C010)+1+EXP(B2*X+C012));
     P1_01=1/(EXP(B0*X+C010)+1+EXP(B2*X+C012));
     P0_11=EXP(B0*X+C110)/(EXP(B0*X+C110)+1+EXP(B2*X+C112));
     P1_11=1/(EXP(B0*X+C110)+1+EXP(B2*X+C112));
     P2_11=EXP(B2*X+C112)/(EXP(B0*X+C110)+1+EXP(B2*X+C112));
     P1_02=1/(EXP(B0*X+C020)+1+EXP(B2*X+C022));
     P1_12=1/(EXP(B0*X+C120)+1+EXP(B2*X+C122));
     P2_12=EXP(B2*X+C122)/(EXP(B0*X+C120)+1+EXP(B2*X+C122));
     P2_22=EXP(B2*X+C222)/(EXP(B0*X+C220)+1+EXP(B2*X+C222));

  OUTPUT;
  END;
RUN;

/*
proc print data=condprob&i;
run;
*/

PROC SORT DATA=CONDPROB&i OUT=ONLYX&i;
BY MTYPE;

ODS TRACE OFF;
ODS OUTPUT FitStatistics=F2stats&i ParameterEstimates=Parm2Ests&i;

PROC LOGISTIC DATA=ONLYX&i(where=&COND) ORDER=DATA;
MODEL MTYPE = &vars2 / LINK=GLOGIT;
WEIGHT W;
RUN;


DATA CONDPARM&i;

  SET Parm2Ests&i END=EOF;
  RETAIN A00 A01 A11 A12 A02 A22 U00 U01 U11 U12 U02 U22;
  IF VARIABLE="X" AND RESPONSE="0" THEN A00=ESTIMATE;
  IF VARIABLE="X" AND RESPONSE="10" THEN A01=ESTIMATE;
  IF VARIABLE="X" AND RESPONSE="11" THEN A11=ESTIMATE;
  IF VARIABLE="X" AND RESPONSE="12" THEN A12=ESTIMATE;
  IF VARIABLE="X" AND RESPONSE="20" THEN A02=ESTIMATE;

  IF VARIABLE="Intercept" AND RESPONSE="0" THEN U00=ESTIMATE;
  IF VARIABLE="Intercept" AND RESPONSE="10" THEN U01=ESTIMATE;
  IF VARIABLE="Intercept" AND RESPONSE="11" THEN U11=ESTIMATE;
  IF VARIABLE="Intercept" AND RESPONSE="12" THEN U12=ESTIMATE;
  IF VARIABLE="Intercept" AND RESPONSE="20" THEN U02=ESTIMATE;
  A22=0;
  U22=0; 
   
  IF EOF THEN OUTPUT;
  DROP ESTIMATE VARIABLE RESPONSE DF STDERR WALDCHISQ PROBCHISQ;
RUN;

data moreparms&i;
  merge condparm&i mparms&i;
run;


DATA UPROBS&i;
  IF _N_=1 THEN SET moreparms&i;
  ELSE SET CONDPROB&i(DROP=&D2) END=EOF;
  IF _N_^=1 THEN DO;

    TOTAL=EXP(A00*X+U00)+EXP(A01*X+U01)+EXP(A11*X+U11)+
          EXP(A12*X+U12)+EXP(A02*X+U02)+EXP(A22*X+U22);

     P00=EXP(A00*X+U00)/TOTAL;
     P01=EXP(A01*X+U01)/TOTAL;
     P11=EXP(A11*X+U11)/TOTAL;
     P12=EXP(A12*X+U12)/TOTAL;
     P02=EXP(A02*X+U02)/TOTAL;
     P22=EXP(A22*X+U22)/TOTAL;

     M2_021=EXP(mom1*x+mom2*x+imprint*x+intmom02)/
               (1+EXP(mom1*x+mom2*x+imprint*x+intmom02));
     M0_021=1/(1+EXP(mom1*x+mom2*x+imprint*x+intmom02));
     M0_010=1/(1+EXP(mom1*x+intmom01));
     M1_010=EXP(mom1*x+intmom01)/(1+EXP(mom1*x+intmom01));
     M0_011=1/(1+EXP(mom1*x+imprint*x+intmom01));
     M1_011=EXP(mom1*x+imprint*x+intmom01)/(1+EXP(mom1*x+imprint*x+intmom01));
     M1_121=1/(1+EXP(mom2*x+imprint*x+intmom12));
     M2_121=EXP(mom2*x+imprint*x+intmom12)/(1+EXP(mom2*x+imprint*x+intmom12));
     M1_122=1/(1+EXP(mom2*x+intmom12));
     M2_122=EXP(mom2*x+intmom12)/(1+EXP(mom2*x+intmom12));



     SUM=P00+P01+P11+P12+P02+P22;

     P000=P00;
     P021=P02*M0_021;
     P201=P02*M2_021;
     P222=P22;
     P010=P01*P0_01*M0_010;
     P011=P01*P1_01*M0_011;
     P100=P01*P0_01*M1_010;
     P101=P01*P1_01*M1_011;
     P110=P11*P0_11;
     P111=P11*P1_11;
     P112=P11*P2_11;
     P121=P12*P1_12*M1_121;
     P122=P12*P2_12*M1_122;
     P211=P12*P1_12*M2_121;
     P212=P12*P2_12*M2_122;
     SUM10=P000+P021+P201+P222+P010+P011+P100+P101+P110+P111+P112+P121+P122+P211+P212;

     OUTPUT;
  END;
RUN;

data combine&i;
  merge mparms&i parms&i;
run;

DATA estimates;
  %IF &i=0 %THEN %STR(SET combine&i;);
  %ELSE %STR(SET estimates combine&i;);
RUN;
*E-STEP;

PROC SORT DATA=UPROBS&i;
BY ID F;

DATA WEIGHTS&i;
  SET UPROBS&i;
  IF ((FLAG_F=0) AND (FLAG_M=0)) THEN OUTPUT;
  IF ((FLAG_F=1) AND (FLAG_M=0)) THEN DO;
     IF ((M=0) AND (C=0) AND (F=0)) THEN DO;
         W=P000/(P000+P010);
         OUTPUT;
         END;
     IF ((M=0) AND (C=0) AND (F=1)) THEN DO;
         W=P010/(P000+P010);
         OUTPUT;
         END;
     IF ((M=0) AND (C=1) AND (F=1)) THEN DO;
         W=P011/(P011+P021);
         OUTPUT;
         END;
     IF ((M=0) AND (C=1) AND (F=2)) THEN DO;
         W=P021/(P011+P021);
         OUTPUT;
         END;
     IF ((M=1) AND (C=0) AND (F=0)) THEN DO;
         W=P100/(P100+P110);
         OUTPUT;
         END;
      IF ((M=1) AND (C=0) AND (F=1)) THEN DO;
         W=P110/(P100+P110);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=0)) THEN DO;
         W=P101/(P101+P111+P121);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=1)) THEN DO;
         W=P111/(P101+P111+P121);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=2)) THEN DO;
         W=P121/(P101+P111+P121);
         OUTPUT;
         END;
     IF ((M=1) AND (C=2) AND (F=1)) THEN DO;
         W=P112/(P112+P122);
         OUTPUT;
         END;
     IF ((M=1) AND (C=2) AND (F=2)) THEN DO;
         W=P122/(P112+P122);
         OUTPUT;
         END;
     IF ((M=2) AND (C=1) AND (F=0)) THEN DO;
         W=P201/(P201+P211);
         OUTPUT;
         END;
     IF ((M=2) AND (C=1) AND (F=1)) THEN DO;
         W=P211/(P201+P211);
         OUTPUT;
         END;
     IF ((M=2) AND (C=2) AND (F=1)) THEN DO;
         W=P212/(P212+P222);
         OUTPUT;
         END;
     IF ((M=2) AND (C=2) AND (F=2)) THEN DO;
         W=P222/(P212+P222);
         OUTPUT;
         END;
  END;
   IF ((FLAG_F=0) AND (FLAG_M=1)) THEN DO;
     IF ((M=0) AND (C=0) AND (F=0)) THEN DO;
         W=P000/(P000+P100);
         OUTPUT;
         END;
     IF ((M=1) AND (C=0) AND (F=0)) THEN DO;
         W=P100/(P000+P100);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=0)) THEN DO;
         W=P101/(P101+P201);
         OUTPUT;
         END;
     IF ((M=2) AND (C=1) AND (F=0)) THEN DO;
         W=P201/(P101+P201);
         OUTPUT;
         END;
     IF ((M=0) AND (C=0) AND (F=1)) THEN DO;
         W=P010/(P010+P110);
         OUTPUT;
         END;
      IF ((M=1) AND (C=0) AND (F=1)) THEN DO;
         W=P110/(P010+P110);
         OUTPUT;
         END;
     IF ((M=0) AND (C=1) AND (F=1)) THEN DO;
         W=P011/(P011+P111+P211);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=1)) THEN DO;
         W=P111/(P011+P111+P211);
         OUTPUT;
         END;
     IF ((M=2) AND (C=1) AND (F=1)) THEN DO;
         W=P211/(P011+P111+P211);
         OUTPUT;
         END;
     IF ((M=1) AND (C=2) AND (F=1)) THEN DO;
         W=P112/(P112+P212);
         OUTPUT;
         END;
     IF ((M=2) AND (C=2) AND (F=1)) THEN DO;
         W=P212/(P112+P212);
         OUTPUT;
         END;
     IF ((M=0) AND (C=1) AND (F=2)) THEN DO;
         W=P021/(P021+P121);
         OUTPUT;
         END;
     IF ((M=1) AND (C=1) AND (F=2)) THEN DO;
         W=P121/(P021+P121);
         OUTPUT;
         END;
     IF ((M=1) AND (C=2) AND (F=2)) THEN DO;
         W=P122/(P122+P222);
         OUTPUT;
         END;
     IF ((M=2) AND (C=2) AND (F=2)) THEN DO;
         W=P222/(P122+P222);
         OUTPUT;
         END;
  END;

DATA iteration&i;
  merge WEIGHTS&i END=EOF;
  BY ID;
  RETAIN LL 0;
  IF ((FLAG_F=0) AND (FLAG_M=0)) THEN DO;
     IF ((M=0) and (F=0) and (C=0)) THEN LL=LL+LOG(P000);
     IF ((M=0) and (F=2) and (C=1)) THEN LL=LL+LOG(P021);
     IF ((M=2) and (F=0) and (C=1)) THEN LL=LL+LOG(P201);
     IF ((M=2) and (F=2) and (C=2)) THEN LL=LL+LOG(P222);
     IF ((M=0) and (F=1) and (C=0)) THEN LL=LL+LOG(P010);
     IF ((M=1) and (F=0) and (C=0)) THEN LL=LL+LOG(P100);
     IF ((M=0) and (F=1) and (C=1)) THEN LL=LL+LOG(P011);
     IF ((M=1) and (F=0) and (C=1)) THEN LL=LL+LOG(P101);
     IF ((M=1) and (F=1) and (C=0)) THEN LL=LL+LOG(P110);
     IF ((M=1) and (F=1) and (C=1)) THEN LL=LL+LOG(P111);
     IF ((M=1) and (F=1) and (C=2)) THEN LL=LL+LOG(P112);
     IF ((M=1) and (F=2) and (C=1)) THEN LL=LL+LOG(P121);
     IF ((M=2) and (F=1) and (C=1)) THEN LL=LL+LOG(P211);
     IF ((M=1) and (F=2) and (C=2)) THEN LL=LL+LOG(P122);
     IF ((M=2) and (F=1) and (C=2)) THEN LL=LL+LOG(P212);
  END;
  IF ((FLAG_F=1) AND (FLAG_M=0) AND LAST.ID) THEN DO;
   IF ((M=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P010);
     IF ((M=0) AND (C=1)) THEN
         LL=LL+LOG(P011+P021);
     IF ((M=1) AND (C=0)) THEN 
         LL=LL+LOG(P100+P110);
     IF ((M=1) AND (C=1)) THEN 
         LL=LL+LOG(P101+P111+P121);         
     IF ((M=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P122);
     IF ((M=2) AND (C=1)) THEN 
         LL=LL+LOG(P211+P201);
     IF ((M=2) AND (C=2)) THEN 
         LL=LL+LOG(P212+P222);  
end;
  IF ((FLAG_F=0) AND (FLAG_M=1) AND LAST.ID) THEN DO;
   IF ((F=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P100);
     IF ((F=0) AND (C=1)) THEN
         LL=LL+LOG(P101+P201);
     IF ((F=1) AND (C=0)) THEN 
         LL=LL+LOG(P010+P110);
     IF ((F=1) AND (C=1)) THEN 
         LL=LL+LOG(P011+P111+P211);         
     IF ((F=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P212);
     IF ((F=2) AND (C=1)) THEN 
         LL=LL+LOG(P121+P021);
     IF ((F=2) AND (C=2)) THEN 
         LL=LL+LOG(P122+P222);  
  END;
  IF EOF THEN 
    OUTPUT;
  KEEP LL;
RUN;

DATA lliters;
  %IF &i=0 %THEN %STR(SET iteration&i;);
  %ELSE %STR(SET lliters iteration&i;);
RUN;


%LET j=&i;

%END;

%MEND;

DATA MISSING;
  SET &YOUR_DATA;
  W=1;
  IF ((FLAG_F=0) AND (FLAG_M=0)) THEN OUTPUT;
  IF FLAG_F=1 THEN DO;
       IF ((M=0) AND (C=0)) THEN DO;
         W=0;
         F=0;
         OUTPUT;
         W=0;
         F=1;
         OUTPUT;
         END;
     IF ((M=0) AND (C=1)) THEN DO;
         W=0;
         F=1;
         OUTPUT;
         W=0;
         F=2;
         OUTPUT;
         END;
     IF ((M=1) AND (C=0)) THEN DO;
         W=0;
         F=0;
         OUTPUT;
         W=0;
         F=1;
         OUTPUT;
         END;
     IF ((M=1) AND (C=1)) THEN DO;
         W=0;
         F=0;
         OUTPUT;
         W=0;
         F=1;
         OUTPUT;
         W=0;
         F=2;
         OUTPUT;
         END;
     IF ((M=1) AND (C=2)) THEN DO;
         W=0;
         F=1;
         OUTPUT;
         W=0;
         F=2;
         OUTPUT;
         END;
     IF ((M=2) AND (C=1)) THEN DO;
         W=0;
         F=0;
         OUTPUT;
         W=0;
         F=1;
         OUTPUT;
         END;
     IF ((M=2) AND (C=2)) THEN DO;
         W=0;
         F=1;
         OUTPUT;
         W=0;
         F=2;
         OUTPUT;
         END;
  END;
  IF FLAG_M=1 THEN DO;
       IF ((F=0) AND (C=0)) THEN DO;
         W=0;
         M=0;
         OUTPUT;
         W=0;
         M=1;
         OUTPUT;
         END;
     IF ((F=0) AND (C=1)) THEN DO;
         W=0;
         M=1;
         OUTPUT;
         W=0;
         M=2;
         OUTPUT;
         END;
     IF ((F=1) AND (C=0)) THEN DO;
         W=0;
         M=0;
         OUTPUT;
         W=0;
         M=1;
         OUTPUT;
         END;
     IF ((F=1) AND (C=1)) THEN DO;
         W=0;
         M=0;
         OUTPUT;
         W=0;
         M=1;
         OUTPUT;
         W=0;
         M=2;
         OUTPUT;
         END;
     IF ((F=1) AND (C=2)) THEN DO;
         W=0;
         M=1;
         OUTPUT;
         W=0;
         M=2;
         OUTPUT;
         END;
     IF ((F=2) AND (C=1)) THEN DO;
         W=0;
         M=0;
         OUTPUT;
         W=0;
         M=1;
         OUTPUT;
         END;
     IF ((F=2) AND (C=2)) THEN DO;
         W=0;
         M=1;
         OUTPUT;
         W=0;
         M=2;
         OUTPUT;
         END;
  END;


/* Creating neccessary indicator variables MTYPE00, MTYPE01, MTYPE11, MTYPE12, MTYPE22, MTYPE02 and variable MTYPE for second polytomous logistic model in equation #2. */

DATA LOGITMOD;
SET MISSING;
  MTYPE01=0;
  MTYPE11=0;
  MTYPE12=0;
  MTYPE00=0;
  MTYPE02=0;
  MTYPE22=0;
  ORDER=0;
  IF (((M=0) AND (F=1)) OR ((M=1) AND (F=0))) THEN DO;
    MTYPE01=1;
    MTYPE=10;
    END;
  IF ((M=1) AND (F=1)) THEN DO;
    MTYPE11=1;
    MTYPE=11;
    END;
  IF (((M=1) AND (F=2)) OR ((M=2) AND (F=1))) THEN DO;
    MTYPE12=1;
    MTYPE=12;
    END;
  IF ((M=0) AND (F=0)) THEN DO;
    MTYPE00=1;
    MTYPE=00;
    END;
  IF (((M=0) AND (F=2)) OR ((M=2) AND (F=0))) THEN DO;
    MTYPE02=1;
    MTYPE=20;
    END;
  IF ((M=2) AND (F=2)) THEN DO;
    MTYPE22=1;
    MTYPE=22;
    END;
    
  IF (M>F) THEN EFFECT=1;
  ELSE IF (M<=F) THEN EFFECT=0;

  IF (((M=0) AND (F=1)) OR ((M=1) AND (F=0))) OR (((M=0) AND (F=2)) OR ((M=2) AND (F=0))) THEN IND1=1;
  ELSE IND1=0;
 
  IF (((M=1) AND (F=2)) OR ((M=2) AND (F=1))) OR (((M=0) AND (F=2)) OR ((M=2) AND (F=0))) THEN IND2=1;
  ELSE IND2=0; 

  INDX1=X*IND1;
  INDX2=X*IND2;
  K1=(C=1);
  K1X=X*(C=1);

  IF (C=1) THEN ORDER=1;
           ELSE ORDER=0;

RUN;

/* Sorting data so that the reference cell is C=1. */

PROC SORT DATA=LOGITMOD OUT=PLOGDATA0;
BY ORDER;
RUN;

%emalg(X MTYPE01 MTYPE11 MTYPE12 MTYPE00 MTYPE02 MTYPE22,X,MTYPE01 MTYPE12 MTYPE02 INDX1 INDX2,k=25);

ods listing;

title2 "Estimates of offspring effects";
proc print data=parm1ests25;
run;

title2 "Estimates of maternal and imprinting effects";
proc print data=mparms1ests25;
run;

proc print data=lliters;
run;
ods listing close;


proc sort data=uprobs25 out=condpr;
by id;
run;

DATA LLWX;
  merge condpr END=EOF;
  BY ID;
  RETAIN LL 0;
  IF ((FLAG_F=0) AND (FLAG_M=0)) THEN DO;
     IF ((M=0) and (F=0) and (C=0)) THEN LL=LL+LOG(P000);
     IF ((M=0) and (F=2) and (C=1)) THEN LL=LL+LOG(P021);
     IF ((M=2) and (F=0) and (C=1)) THEN LL=LL+LOG(P201);
     IF ((M=2) and (F=2) and (C=2)) THEN LL=LL+LOG(P222);
     IF ((M=0) and (F=1) and (C=0)) THEN LL=LL+LOG(P010);
     IF ((M=1) and (F=0) and (C=0)) THEN LL=LL+LOG(P100);
     IF ((M=0) and (F=1) and (C=1)) THEN LL=LL+LOG(P011);
     IF ((M=1) and (F=0) and (C=1)) THEN LL=LL+LOG(P101);
     IF ((M=1) and (F=1) and (C=0)) THEN LL=LL+LOG(P110);
     IF ((M=1) and (F=1) and (C=1)) THEN LL=LL+LOG(P111);
     IF ((M=1) and (F=1) and (C=2)) THEN LL=LL+LOG(P112);
     IF ((M=1) and (F=2) and (C=1)) THEN LL=LL+LOG(P121);
     IF ((M=2) and (F=1) and (C=1)) THEN LL=LL+LOG(P211);
     IF ((M=1) and (F=2) and (C=2)) THEN LL=LL+LOG(P122);
     IF ((M=2) and (F=1) and (C=2)) THEN LL=LL+LOG(P212);
  END;
   IF ((FLAG_F=1) AND (FLAG_M=0) AND LAST.ID) THEN DO;
   IF ((M=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P010);
     IF ((M=0) AND (C=1)) THEN
         LL=LL+LOG(P011+P021);
     IF ((M=1) AND (C=0)) THEN 
         LL=LL+LOG(P100+P110);
     IF ((M=1) AND (C=1)) THEN 
         LL=LL+LOG(P101+P111+P121);         
     IF ((M=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P122);
     IF ((M=2) AND (C=1)) THEN 
         LL=LL+LOG(P211+P201);
     IF ((M=2) AND (C=2)) THEN 
         LL=LL+LOG(P212+P222);  
  END;  
   IF ((FLAG_F=0) AND (FLAG_M=1) AND LAST.ID) THEN DO;
   IF ((F=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P100);
     IF ((F=0) AND (C=1)) THEN
         LL=LL+LOG(P101+P201);
     IF ((F=1) AND (C=0)) THEN 
         LL=LL+LOG(P010+P110);
     IF ((F=1) AND (C=1)) THEN 
         LL=LL+LOG(P011+P111+P211);         
     IF ((F=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P212);
     IF ((F=2) AND (C=1)) THEN 
         LL=LL+LOG(P121+P021);
     IF ((F=2) AND (C=2)) THEN 
         LL=LL+LOG(P122+P222);  
  END;
 IF EOF THEN 
    OUTPUT;
  KEEP LL;
RUN;


%emalg(X MTYPE01 MTYPE11 MTYPE12 MTYPE00 MTYPE02 MTYPE22,X,MTYPE01 MTYPE12 MTYPE02,k=25);

ods listing;

proc print data=lliters;
run;
ods listing close;

proc sort data=uprobs25 out=condprwox;
by id;
run;

DATA LLWOX;
  merge condprwox END=EOF;
  BY ID;
  RETAIN LL 0;
   IF ((FLAG_F=0) AND (FLAG_M=0)) THEN DO;
     IF ((M=0) and (F=0) and (C=0)) THEN LL=LL+LOG(P000);
     IF ((M=0) and (F=2) and (C=1)) THEN LL=LL+LOG(P021);
     IF ((M=2) and (F=0) and (C=1)) THEN LL=LL+LOG(P201);
     IF ((M=2) and (F=2) and (C=2)) THEN LL=LL+LOG(P222);
     IF ((M=0) and (F=1) and (C=0)) THEN LL=LL+LOG(P010);
     IF ((M=1) and (F=0) and (C=0)) THEN LL=LL+LOG(P100);
     IF ((M=0) and (F=1) and (C=1)) THEN LL=LL+LOG(P011);
     IF ((M=1) and (F=0) and (C=1)) THEN LL=LL+LOG(P101);
     IF ((M=1) and (F=1) and (C=0)) THEN LL=LL+LOG(P110);
     IF ((M=1) and (F=1) and (C=1)) THEN LL=LL+LOG(P111);
     IF ((M=1) and (F=1) and (C=2)) THEN LL=LL+LOG(P112);
     IF ((M=1) and (F=2) and (C=1)) THEN LL=LL+LOG(P121);
     IF ((M=2) and (F=1) and (C=1)) THEN LL=LL+LOG(P211);
     IF ((M=1) and (F=2) and (C=2)) THEN LL=LL+LOG(P122);
     IF ((M=2) and (F=1) and (C=2)) THEN LL=LL+LOG(P212);
  END;
  IF ((FLAG_F=1) AND (FLAG_M=0) AND LAST.ID) THEN DO;
   IF ((M=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P010);
     IF ((M=0) AND (C=1)) THEN
         LL=LL+LOG(P011+P021);
     IF ((M=1) AND (C=0)) THEN 
         LL=LL+LOG(P100+P110);
     IF ((M=1) AND (C=1)) THEN 
         LL=LL+LOG(P101+P111+P121);         
     IF ((M=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P122);
     IF ((M=2) AND (C=1)) THEN 
         LL=LL+LOG(P211+P201);
     IF ((M=2) AND (C=2)) THEN 
         LL=LL+LOG(P212+P222);
  END;
  IF ((FLAG_F=0) AND (FLAG_M=1) AND LAST.ID) THEN DO;
   IF ((F=0) AND (C=0)) THEN 
         LL=LL+LOG(P000+P100);
     IF ((F=0) AND (C=1)) THEN
         LL=LL+LOG(P101+P201);
     IF ((F=1) AND (C=0)) THEN 
         LL=LL+LOG(P010+P110);
     IF ((F=1) AND (C=1)) THEN 
         LL=LL+LOG(P011+P111+P211);         
     IF ((F=1) AND (C=2)) THEN 
         LL=LL+LOG(P112+P212);
     IF ((F=2) AND (C=1)) THEN 
         LL=LL+LOG(P121+P021);
     IF ((F=2) AND (C=2)) THEN 
         LL=LL+LOG(P122+P222);  
  END;
IF EOF THEN  OUTPUT;
  KEEP LL;
RUN;

DATA EMMISSLR;
  MERGE llwx(RENAME=(LL=WX)) 
        llwox(RENAME=(LL=WOX));
  LRT=-2*(WOX-WX);
  ProbChiSq=1-PROBCHI(LRT,2);
RUN;

ODS LISTING;
OPTIONS NOTES;

proc print data=emmisslr;
run;