library(boot) library(MASS) library(base) library(stats) ######################################################################### # Score Test and Improved Score Test. y(vector) is trait value, # # and each row of x represents the genotypic score of one individual. # ######################################################################### imp_score_test=function(y,x) { e=10^(-6) N=length(y) y=y-mean(y) a=sum(y*y)/N if(N==length(x)) { x=x-mean(x) M=1 U=sum(x*y) V1=sum(x*x) t1=U^2/(a*V1+e) t2=U^2/(a*V1-U^2/N) } else { M=ncol(x) Y=matrix(y,N,1) x=t(t(x)-colMeans(x)) U=t(x)%*%Y t1=t(U)%*%ginv(t(x)%*%x)%*%U/a t2=t(U)%*%ginv(t(x)%*%x*a-U%*%t(U)/N)%*%U } p1=1-pchisq(t1,M) p2=1-pchisq(t2,M) pvalue=c(p1,p2) } ######################################################################################### # Score Test and Improved Score Test that are robust to population # # stratification. y(vector) repesents trait value. Each row of x represents the # # genotypic score of one individual. Each row of G represents the genotypic score # # of one individual at a set of unlinked markers. # ######################################################################################### imp_score_testsp=function(y,x,G) { K=10 N=length(y) L=ncol(G) G=t(t(G)-colMeans(G)) if(N>=L) { A=t(G)%*%G B=eigen(A) E=B$vectors T=G%*%E[,1:K] } else { GG=G%*%t(G) BB=eigen(GG) EE=BB$vectors T=EE[,1:K] } d=diag(t(T)%*%T) Beta=t(T)%*%y/d y1=y-mean(y)-T%*%Beta x1=x if(length(x)==N) { Beta=t(T)%*%x/d x1=x-mean(x)-T%*%Beta } else { M=ncol(x) for(i in 1:M) { Beta=t(T)%*%x[,i]/d x1[,i]=x[,i]-mean(x[,i])-T%*%Beta } } pvalue=imp_score_test(y1,x1) }