Based on chapter 12.3.2 and implementing Higham's algorithm, I simulated the following:
/******************
* test simulation *
******************/
* step 1: create 5 random SNPS;
**********************************;
* in practice you retrieve 5 snps that are located nearby and you assemble them in a dataframe;
data Categories(keep=V1--V5);
call streaminit(4321);
array p[3] (0.4 0.2 0.4); /* probabilities */
do i = 1 to 100;
V1 = rand("Table", of p[*]); /* use OF operator */
V2 = rand("Table", of p[*]); /* use OF operator */
V3 = rand("Table", of p[*]); /* use OF operator */
V4 = rand("Table", of p[*]); /* use OF operator */
V5 = rand("Table", of p[*]); /* use OF operator */
output;
end;
run;
proc print data=Categories;
run;
data snps (keep= snp1--snp5);
set Categories;
snp1 = V1-1;
snp2 = V2-1;
snp3 = V3-1;
snp4 = V4-1;
snp5 = V5-1;
keep snp1--snp
run;
proc print data=snps;
run;
proc freq data=snps;
tables snp1 snp2 snp3 snp4 snp5/ nocum nopercent;
run;
* step2: create a kinship matrix, I will use the AR1 structure to get a quick sort of similarity matrix;
*******************************************************************************************;
proc iml;
/* return p x p matrix whose (i,j)th element is rho^|i - j| */
start AR1Corr(rho, n);
return rho##distance(T(1:n), T(1:n), "L1");
finish;
/* test on 100 x 100 matrix with rho = 0.8 */
rho = 0.8;
n = 100;
kinship = AR1Corr(rho, n);
print kinship[format=Best7.];
*export to dataset;
create dataK from kinship;
append from kinship;
close dataK;
quit;
* simulate a trait;
*******************;
* in practice read the kinship from csv file, create matrix from dataset;
proc iml;
* import kinship from data;
use work.dataK; /* open the data set */
read all var _ALL_ into K;
close work.dataK;
print(K);
* implement Higham's algorithm for computing the nearest correlation matrix to a given symmetric matrix;
* Ref: https://blogs.sas.com/content/iml/2012/11/28/computing-the-nearest-correlation-matrix.html ;
/* Project symmetric X onto S={positive semidefinite matrices}.
Replace any negative eigenvalues of X with zero */
start ProjS(X);
call eigen(D, Q, X); /* note X = Q*D*Q` */
V = choose(D>0, D, 0);
W = Q#sqrt(V`); /* form Q*diag(V)*Q` */
return( W*W` ); /* W*W` = Q*diag(V)*Q` */
finish;
/* project square X onto U={matrices with unit diagonal}.
Return X with the diagonal elements replaced by ones. */
start ProjU(X);
n = nrow(X);
Y = X;
diagIdx = do(1, n*n, n+1);
Y[diagIdx] = 1; /* set diagonal elements to 1 */
return ( Y );
finish;
/* Helper function: the infinity norm is the max abs value of the row sums */
start MatInfNorm(A);
return( max(abs(A[,+])) );
finish;
/* Given a symmetric correlation matrix, A,
project A onto the space of positive semidefinite matrices.
The function uses the algorithm of Higham (2002) to return
the matrix X that is closest to A in the Frobenius norm. */
start NearestCorr(A);
maxIter = 100; tol = 1e-8; /* parameters...you can change these */
iter = 1; maxd = 1; /* initial values */
Yold = A; Xold = A; dS = 0;
do while( (iter <= maxIter) & (maxd > tol) );
R = Yold - dS; /* dS is Dykstra's correction */
X = ProjS(R); /* project onto S={positive semidefinite} */
dS = X - R;
Y = ProjU(X); /* project onto U={matrices with unit diagonal} */
/* How much has X changed? (Eqn 4.1) */
dx = MatInfNorm(X-Xold) / MatInfNorm(X);
dy = MatInfNorm(Y-Yold) / MatInfNorm(Y);
dxy = MatInfNorm(Y - X) / MatInfNorm(Y);
maxd = max(dx,dy,dxy);
iter = iter + 1;
Xold = X; Yold = Y; /* update matrices */
end;
return( X ); /* X is positive semidefinite */
finish;
KK = NearestCorr(K);
/* for large matrices, might need to correct for numerical rounding errors */
epsilon = 1e-10;
KKK = ProjU( KK/(1+epsilon) ); /* divide off-diagonal elements by 1+epsilon */
* import snps from data;
use work.snps; /* open the data set */
read all var _ALL_ into GT;
close work.snps;
print(GT);
snp1 = GT[,1];
snp2 = GT[,1];
snp3 = GT[,1];
snp4 = GT[,1];
snp5 = GT[,1];
* set parms;
n=100; /* 100 subjects */
sigma2_g = 0.1; /* Parameter 1: residual covariance */
sigma2_e = 0.5;
beta = {1, 0.5, 0.7, 3, 0.7, 0.5}; /* ratio data, value of 1 for intercept*/
* varcov;
G = sigma2_g*KKK ;
R = sigma2_e*I(n);
varcov = G + R ;
* subjectIDs;
treeID = colvec(1:n);
*print(treeID);
* design matrix;
X = J(n,1,1)||GT ; /*col vector of 1 for intercept*/
* simulate trait;
call randseed(1234);
create Mix var {"treeID" "Trait1" "snp1" "snp2" "snp3" "snp4" "snp5"};
/* create dataset testsim */
* Note: if you are trying to simulate random multivariate normal data,
you must use a positive definite matrix.;
zero = j(1, n, 0); /* the zero vector */
eps = RandNormal(1, zero, varcov); /* eps ~ MVN(0,R) */
trait1 = X*beta + eps`; /* fixed effects, correlated errors */
append;
close;
quit;
* analysis;
**************;
proc sgplot data=Mix;
title "distribution parameters";
histogram trait1;
density trait1;
density trait1/type=kernel;
run;
* you clearly see a biodal distribution;
* analysis;
**********************************;
* the kinship data set must have in
- column 1 : the treeIDs
- column 2 : parm (all 1's)
- column 3 : values from 1 to n ;
data kin ;
retain treeID parm row ;
set datak;
treeID = _N_ ;
parm = 1 ;
row = _N_ ;
run;
Proc mixed data=Mix METHOD=REML noprofile covtest itdetails ;
class treeID ;
model trait1 = snp1 snp2 snp3 snp4 snp5 /solution ddfm=KR chisq;
random treeID / type=lin(1) ldata=work.kin ;
parms/ lowerb= 0 0;
ods output covParms=myParms;
run;
I do however not have the expected results. I can't see what I did wrong.
Despite the use of Higham's algorithm, I get the notes:
An infinite likelihood is assumed in iteration 1 because of a nonpositive definite estimated V matrix for Subject 1. NOTE: Convergence criteria met. NOTE: Estimated G matrix is not positive definite.
... View more
