Sample SAS for Real Data Set

  /*================================================================*
 *This SAS program is to analyze twin data. It has three sections.*
 *The first section reads the data. You need to specify your own  *
 *path for the data set and the variable names to be input. This  *
 *is also where you create your new variables.                    *
 *The second section should not be changed. This section creates  *
 *random effects to accommodate different correlations of MZ and  *
 *DZ twins.                                                       *
 *The third section executes different models. You need to include*
 *appropriate variables for your analysis as well as try different*
 *mode of inheritance and initial parameter values.               *
 *Please email heping.zhang@yale.edu if you discovered any errors *
 *with this program. This program can be freely used and          *
 *distributed. However, the owners of this program is not         *
 *responsible for inaccuracy or errors as a result of using this  *
 *program. It is the responsibility of the users to verify the    *
 *the validity of this program and conduct their analysis         * 
 *appropriately.                                                  *
 *================================================================*/

/*Begin Section 1*/
data one;
*you need to enter your own file path in the next statement;
INFILE 'D:\shared\hz2\Students\Gongfu\biometrics-data.txt' missover delimiter='09'x  firstobs=2;
*you need to include appropriate variables the next statement;
input FAMID     Sex     GA      BW      RDS     BPD     INST zygo;
*you need to create appropriate additional variables if needed;
inst2=(INST=2);
inst3=(INST=3);
id=_n_;
run;
/*End Section 1*/
/*Begin Section 2: no changes are needed in this section*/
proc iml;
covA1=sqrsym({1,1,1});
gA1=t(root(covA1));
covA2=sqrsym({1,0.5,1});
gA2=t(root(covA2));
covD1=sqrsym({1,1,1});
gD1=t(root(covD1));
covD2=sqrsym({1,0.25,1});
gD2=t(root(covD2));

/* for additive genetic effect, the G design matrix determines the model through
the model covariance structure: VAR[A]=GZG`. The mixed model is then represented
as: logit(p{y=1}) = mu + ZG + e */

g1=INSERT(gA1,gD1,0,3);                 /*combine matrix for MZ twin*/
g2=INSERT(gA2,gD2,0,3);                 /*combine matrix for DZ twin*/
use one;
        read all var {FAMID zygo} into vzygo;
close one;
nsub=nrow(vzygo);
g=j(nsub,4,.);
do i=1 to nsub;
        ind=2-mod(i,2);
        if(vzygo[i,2]=1) then g[i,]=g1[ind,];                   /* MZ twin*/
        if(vzygo[i,2]=0) then g[i,]=g2[ind,];                   /* DZ twin*/
end;
cname={"aone" "atwo" "done" "dtwo"};
create rmatrix from g[colname=cname];
append from g;
close rmatrix;
quit;
data rmatrix;
set rmatrix;
id=_n_;
run;
data two;
merge one rmatrix;
by id;
run;
/*End Section 2*/
/*Begin Section 3*/
title1 'ADE model';
proc nlmixed data=two;
*you need to choose your own initial values;
*it may be useful to run a simple model such as PROC LOGISTIC;
*usually you may need to try-and-err;
parms beta1=-0.192 beta2=0.088 beta3=0.0011 beta4=0.889 beta5=0.230  beta6=-1.092 s1=4 s2=0.01;
*include the fixed covariates in the next statement;
fixed1 = beta0+beta1*SEX+beta2*GA+beta3*BW+beta4*inst2+beta5*inst3+beta6*RDS;
*the following statement specifies additive (a1 and a2) and;
*dominant (d1 and d2) random effects;
random1 = a1*aone + a2*atwo + d1*done + d2*dtwo;
eta = fixed1 + random1;
*specifies a probit model;
p = 1 - probnorm(eta);
model BPD ~ binary(p);
random a1 a2 d1 d2 ~ normal([0,0,0,0],[s1,0,s1,0,0,s2,0,0,0,s2]) subject=FAMID;
estimate 'icc' (s1+s2)/(1+s1+s2);
run;

title2 'ACE model';
proc nlmixed data=two;
parms beta0=2.793 beta1=0.133 beta2=-0.095 beta3=-0.0011 beta4=-0.929 beta5=-0.283  beta6=1.112 s1=3 s2=1;
fixed1 = beta0+beta1*SEX+beta2*GA+beta3*BW+beta4*inst2+beta5*inst3+beta6*RDS;
*the following statement specifies additive (a1 and a2) and;
*common environment (c) random effects;
random1 = a1*aone+a2*atwo+c;
eta = fixed1 + random1;
p = 1 - probnorm(eta);
model BPD ~ binary(p);
random a1 a2 c ~ normal([0,0,0],[s1,0,s1,0,0,s2]) subject=FAMID;
estimate 'icc' s1/(1+s1+s2);
run;
/*End Section 3*/