######################################################## ## MHT-O is an association test to test association ## ## between multiple traits and a genetic variant. x ## ## represents genotype (number of minor allele) ## ## vector. y represents trait values (a matrix). ## ## Each row of y represents an individual and each ## ## column of y represents a trait. numper represents ## ## the number of permutations. ## ######################################################## MHT_O=function(x,y,numper) { n=ncol(y) N=nrow(y) a=lm(y~x) aa=sqrt(var(x)) u=(a$coeff[2,])*aa res=a$resid varr=var(res) T=ginv(u%*%t(u)+varr) w=rep(0,n) w[n]=as.numeric(t(u)%*%T%*%u) for (i in 1:(n-1)) { #writeLines(paste("w=",w,"i=",i)) if (i==1) a=Max(u,res) # else a=Max(a[[2]],a[[3]]) w[n-i]=a[[1]] # else break # writeLines(paste("a=",a[[1]])) } ###################### wper=matrix(0,numper,n) for(j in 1:numper) { x=sample(x) a=lm(y~x) u=(a$coeff[2,])*aa res=a$resid varr=var(res) T=ginv(u%*%t(u)+varr) wper[j,n]=as.numeric(t(u)%*%T%*%u) for (i in 1:(n-1)) { if (i==1) a=Max(u,res) else a=Max(a[[2]],a[[3]]) wper[j,n-i]=a[[1]] #else break } } #pv=sum(wper>w)/numper+sum(wper==w)/numper/2 T=rbind(w,wper) T0=apply(T,2,rank) T=apply(T0,1,max) pv=sum(T[-1]>T[1])/numper+sum(T[-1]==T[1])/numper/2 pv } ############################################### Max=function(u,res) { n=length(u) w=rep(0,n) for (i in 1:n) { varr=var(res[,-i]) T=ginv(u[-i]%*%t(u[-i])+varr) w[i]=as.numeric(t(u[-i])%*%T%*%u[-i]) } ind=which(w==max(w))[1] list(max(w),u[-ind],res[,-ind]) } ####################################################