library(boot) library(MASS) library(base) library(stats) library(Matrix) ################################################################## ## MultP-PE is a statistical method to test associations ## ## between a genetic variant and Multiple Phenotypes ## ## based on cross-validation Prediction Error. ## ## ## ## In MultP-PE function, ## ## x: represents genotype(number of minor alleles) matrix. ## ## Each row represents an individual and each column ## ## represents a SNP. ## ## y: represents phenotype matrix. Each row represents ## ## an individual and each column represents a phenotype. ## ## numper: represents the number of permutations to ## ## obtain the p-value of the test statistic. ## ## lambda: represents the penalized parameters in ridge ## ## regression(vector). ## ################################################################## MultP-PE=function(x,y,numper,lambda) { n=nrow(y) p=ncol(x) yf=cbind(rep(1,n),y) mf=ncol(yf) nlambda=length(lambda) J=matrix(1,n,nlambda) SVD=svd(yf,nu=ncol(yf),nv=ncol(yf)) U=SVD$u V=SVD$v d=SVD$d H=matrix(0,n,nlambda) C=matrix(0,mf,nlambda) pv=matrix(0,p,1) for(pp in 1:p) { xx=x[,pp] xm=t(U)%*%xx for(i in 1:nlambda) { C[,i]=as.matrix(d^2/(d^2+lambda[i])) H[,i]=rowSums(U*t(t(U)*as.vector(C[,i]))) } xhat_lambda=U%*%(C*as.vector(xm)) B=(as.vector(xx)-xhat_lambda)/(J-H) T0=-colSums(B*B) ############################# Tper=matrix(0,numper,nlambda) for(j in 1:numper) { ind=sample(n) x1=xx[ind] xm=t(U)%*%x1 xhat_lambda=U%*%(C*as.vector(xm)) B=(as.vector(x1)-xhat_lambda)/(J-H) Tper[j,]=-colSums(B*B) } TT=rbind(T0,Tper) Tr=apply(TT,2,rank) Tv=apply(Tr,1,max) pv[pp]=sum(Tv[-1]>Tv[1])/numper+sum(Tv[-1]==Tv[1])/numper/2 } pv }