/* READ COMMENTS:
This is the missing data additive logistic approach for multiple offspring with missing fathers or mothers.  

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


This code allows for up to three children per family.  

The file lib2.iml must also be downloaded.  You must change the line reading 

%INCLUDE "..\IML\LIB2.IML"/ SOURCE2;

to point to the location where lib2.iml is saved.  

Here 3 imputations are computed.  
It is recommended that more iterations are tried to confirm the results if more than 30% of the parents are missing. 

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 ith child carries, and
X= the quantitative trait for the ith child.
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.  
FAMILYID= family identifier.
CHILDID= identifier of which child in the family. 
NKIDS= number of children in the family.

YOU MUST ENTER THE NAME OF YOUR DATA SET IN THE FIRST %LET STATEMENT BELOW!

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


TITLE1 "POLYTOMOUS logistic approach for missing fathers or mothers and multiple offspring";

%LET YOUR_DATA=;

ODS LISTING CLOSE;

DATA MATINGTYPES;
  SET &YOUR_DATA;
  MTYPE01=0;
  MTYPE11=0;
  MTYPE12=0;
  MTYPE00=0;
  MTYPE02=0;
  MTYPE22=0;
  IF (((M=0) AND (F=1)) OR ((M=1) AND (F=0))) THEN DO;
    MTYPE01=1;
    FAMILY=1112;
    END;
  IF ((M=1) AND (F=1)) THEN DO;
    MTYPE11=1;
    FAMILY=1212;
    END;
  IF (((M=1) AND (F=2)) OR ((M=2) AND (F=1))) THEN DO;
    MTYPE12=1;
    FAMILY=1222;
    END;
  IF ((M=0) AND (F=0)) THEN DO;
    MTYPE00=1;
    FAMILY=1111;
    END;
  IF (((M=0) AND (F=2)) OR ((M=2) AND (F=0))) THEN DO;
    MTYPE02=1;
    FAMILY=1122;
    END;
  IF ((M=2) AND (F=2)) THEN DO;
    MTYPE22=1;
    FAMILY=2222;
    END;
  if ((flag_f=1) or (flag_m=1)) then missing=1;
  else missing=0;
;

proc means data=matingtypes;
var nkids;
output out=mostkids max=maxkids;


data maxvar;
  SET mostkids(drop=_type_ _freq_);
  call symput('max', trim(left(maxkids)));
run;

%PUT &max;



data array;
  IF _N_=1 THEN SET mostkids(drop=_type_ _freq_);
  ELSE SET matingtypes END=EOF;
  by familyid;
  IF _N_^=1 then do;
  retain k1-k&max qt1-qt&max;
  array kids{&max} k1-k&max; 
  array qts{&max} qt1-qt&max;
  if first.familyid then do;
    do j=1 to &max;
       kids{j}=.;
       qts{j}=.;
    end;
  end;
  do j=1 to &max;
    if j=childid then do;
       kids{j}=c;
       qts{j}=x;
       if last.familyid then output;
       drop j x c i famadd childid order;
    end;
  end;
  END;
run;

data dropmiss;
  set array;
  leaveout=0;
  if (flag_f=1) then do;
    if (nkids=1) then do;
      if (((m=0) and (k1=0)) or 
          ((m=2) and (k1=2))) then leaveout=1;
    end;
    if (nkids=2) then do;
      if (((m=0) and (k1=0) and (k2=0)) or
          ((m=2) and (k1=2) and (k2=2))) then leaveout=1;
    end;
    if (nkids=3) then do;
      if (((m=0) and (k1=0) and (k2=0) and (k3=0)) or
          ((m=2) and (k1=2) and (k2=2) and (k3=2))) then leaveout=1;
    end;
  end;
  if (flag_m=1) then do;
    if (nkids=1) then do;
      if (((f=0) and (k1=0)) or 
          ((f=2) and (k1=2))) then leaveout=1;
    end;
    if (nkids=2) then do;
      if (((f=0) and (k1=0) and (k2=0)) or
          ((f=2) and (k1=2) and (k2=2))) then leaveout=1;
    end;
    if (nkids=3) then do;
      if (((f=0) and (k1=0) and (k2=0) and (k3=0)) or
          ((f=2) and (k1=2) and (k2=2) and (k3=2))) then leaveout=1;
    end;
  end;
  if leaveout=0 then output dropmiss;

run;


proc iml;

%INCLUDE "..\IML\LIB2.IML"/ SOURCE2;

USE dropmiss; 
if (&max=1) then do;                 
READ ALL VAR {m f k1 qt1 nkids family missing flag_f flag_m familyid};   
end;
if (&max=2) then do;                 
READ ALL VAR {m f k1 k2 qt1 qt2 nkids family missing flag_f flag_m familyid};   
end;
if (&max=3) then do;                 
READ ALL VAR {m f k1 k2 k3 qt1 qt2 qt3 nkids family missing flag_f flag_m familyid};   
end;
  /* read in type of family and associated vector of cumulative sums */
CLOSE dropmiss;

nparms=14;
nparms2=15;
additive=4;
totalbeta=j(additive,1,0);
totalcovs=j(additive,additive,0);
qpldf=2;
adddf=1;

if (&max=1) then do;
x=qt1;
kid=k1;
end;
if (&max=2) then do;
x=(qt1||qt2);
kid=(k1||k2);
end;
if (&max=3) then do;
x=(qt1||qt2||qt3);
kid=(k1||k2||k3);
end;

ods listing;



run add_will(x,kid,family,nkids,missing,additive,addlambda0,addcovwill);
run mpl_will(x,kid,family,nkids,missing,nparms2,mlambda0);



run addmiss(x,kid,m,f,family,flag_m,flag_f,nkids,missing,addlambda0,mlambda0,missm,missf,mfamily,notmissing);

run add_will(x,kid,mfamily,nkids,notmissing,additive,addlambda1,addcovwill);
totalbeta=totalbeta+addlambda1;
totalcovs=totalcovs+addcovwill;



run addmiss(x,kid,m,f,family,flag_m,flag_f,nkids,missing,addlambda0,mlambda0,missm,missf,mfamily,notmissing);

run add_will(x,kid,mfamily,nkids,notmissing,additive,addlambda2,addcovwill);
totalbeta=totalbeta+addlambda2;
totalcovs=totalcovs+addcovwill;



run addmiss(x,kid,m,f,family,flag_m,flag_f,nkids,missing,addlambda0,mlambda0,missm,missf,mfamily,notmissing);

run add_will(x,kid,mfamily,nkids,notmissing,additive,addlambda3,addcovwill);
totalbeta=totalbeta+addlambda3;
totalcovs=totalcovs+addcovwill;



avebeta3=totalbeta/3;
avecov3=totalcovs/3;
estcov=j(additive,additive,0);
halfvar=addlambda1-avebeta3;
vars=halfvar*halfvar`;
estcov=estcov+vars;
halfvar=addlambda2-avebeta3;
vars=halfvar*halfvar`;
estcov=estcov+vars;
halfvar=addlambda3-avebeta3;
vars=halfvar*halfvar`;
estcov=estcov+vars;
estB=estcov*(2/3);
varcov=avecov3+estB;
test=j(adddf,additive,0);
test[1,1]=1;
halftest=test*varcov;
testcov=halftest*test`;
covinv=inv(testcov);
lavebeta=test*avebeta3;
halfwald=covinv*lavebeta;
waldtest=lavebeta`*halfwald;
pvalue=1-probchi(waldtest,adddf);
ods listing;
print waldtest;
ods listing close;