# Sample SAS Program

` ````
/*================================================================*
*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 'Z:\yj56\twin studies\080312P\biometrics-simb.txt' missover delimiter='09'x firstobs=2;
*you need to include appropriate variables the next statement;
input FAMID BPD SEX X1 ZYGO ID;
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=gA1||gD1; /*combine matrix for MZ twin*/
g2=gA2||gD2; /*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.1, beta2=-0.79, s1=0.1, s3=0.1;
*include the fixed covariates in the next statement;
fixed1 = beta0+beta1*SEX+beta2*X1;
*the following statement specifies additive (a1 and a2) and;
*common shared (c) random effects;
random1 = a1*aone+a2*atwo+c;
eta = fixed1 + random1;
*specifies a probit model;
p = probnorm(eta);
model BPD ~ binary(p);
random a1 a2 c ~ normal([0,0,0],[s1,0,s1,0,0,s3]) subject=FAMID;
estimate 'heritability' s1/(1+s1+s3);
run;
title2 'ACE model';
proc nlmixed data=two;
parms beta1 = 0.1, beta2=-0.79, s1=0.1, s2=0.1;;
fixed1 = beta0+beta1*SEX+beta2*X1;
*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;
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 'heritability' (s1+s2)/(1+s1+s2);
run;
/*End Section 3*/
```